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ABSTRACT 

The  problem  of  narrow  band  interference  while  transmitting  broad-band  signals 
like  Direct  Sequence  Spread  Spectrum  is  a  common  source  of  problems  in  Electronic 
Warfare.  This  can  occur  either  due  to  intentional  jamming  or  due  to  unavoidable 
signal  sources  present  in  the  vicinity  of  the  receiver.  Lack  of  improper  information 
on  these  narrow  band  interferers  makes  it  difficult  to  cancel  them. 

In  this  thesis  the  above  problem  is  addressed  by  using  an  adaptive  notch  filtering 
technique.  Before  adopting  such  a  technique  other  methods  like  the  Least  Square 
Estimator  and  the  Maximum  Likelihood  Estimator  were  explored.  However  the  Kwan 
and  Martin  adaptive  notch  filter  structure  was  found  both  relevant  and  suitable  for 
the  problem  of  interest.  The  Kwan  and  Martin  method  has  the  difficulty  of  increasing 
hardware  complexity  with  number  of  notches.  This  makes  it  difficult  to  implement 
in  real  time.  A  new  algorithm  was  developed  for  the  purpose  of  implementing  the 
structure  in  real-time.  This  new  algorithm  offers  the  same  performance  at  reduced 
hardware  complexity.  This  algorithm  was  simulated  and  the  results  were  presented.  A 
hardware  feasibility  is  discussed  by  proposing  a  simple  structure  based  upon  existing 
commercial  signal  processing  chips. 
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I.  INTRODUCTION 

A.      MOTIVATION 

The  problem  of  suppressing  narrow-band  noise  from  a  broad-band  spectrum  is 
of  considerable  importance  in  the  areas  of  Electronic  Warfare  and  Anti  Submarine 
Warfare  scenarios.  In  this  thesis  a  workable  solution  is  proposed  for  the  problem  cre- 
ated by  narrow-band  interference.  This  solution  will  come  in  the  form  of  an  adaptive 
digital  notch  filter.  However  before  the  details  of  the  solution  to  the  problem  are 
discussed,  a  brief  description  of  the  problem  is  given 

1.  Electronic  Warfare(EW) 

A  typical  EW  scenario  is  given  in  Fig  1.1.  In  this  scenario  a  transmitter 
transmits  information  over  a  long  distance.  As  an  Electronic  Counter  Measure(ECM) 
this  information  is  coded  and  transmitted  as  a  Direct  Sequence  Spread  Spectrum 
(DSSS)  signal.  Inherently  this  signal  is  a  broad-band  signal.  However  at  the  EW 
receiver  there  is  often  narrow-band  interference  from  such  things  as  push  to  talk 
systems(PTS)  that  swamp  out  the  received  DSSS  signal.  Sometimes  for  reasons  of 
signal  security  PTS  operating  frequencies  are  varying  with  time.  Under  conditions 
such  as  these  the  EW  receiver  cannot  function  effectively.  For  proper  functioning  of 
the  EW  receiver  we  must  enhance  the  received  signal  by  selectively  and  adaptively 
suppressing  these  narrow-band  signals. 

2.  Assumptions 

This  thesis  addresses  the  problem  arising  at  the  tactical  data  communica- 
tion link.  In  this  EW  scenario  as  dipicted  in  Fig  1.1  we  are  attempting  to  perform 
signal  analysis  on  a  DSSS  signal  using  a  EW  reciever.  This  EW  receiver  is  not  the  des- 
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ignated  receiver  and  hence  does  not  have  the  spreading  code  available.  For  this  reason 
the  narrow-band  interference  severely  limits  the  ability  to  determine  the  charaterstics 
of  the  recieved  signal  such  as  carrier  frequency,  chip-frequency  etc. 

A  typical  DSSS  system  is  the  Tactical  Information  Distribution  Systems 
(JTIDS)[Ref.  7].  Data  on  such  a  link  is  generally  in  the  form  of  digital  messages. 
A  typical  L-band  (960-1215  Mhz)  JTIDS  allows  transfer  of  digital  data  between  any 
properly  equipped  users  within  line-of-sight.  The  typical  RF  band-width  required  is 
around  lOMhz. 

Assuming  a  Superheterodyne  EW-receiver,  the  band-width  of  the  Interme- 
diate Frequency  (IF)  need  be  only  lOMhz.  This  enables  us  to  digitize  the  signal  at 
the  base-band  with  a  (20-50)Mhz  Analog  to  Digital(A/D)  converter.  It  is  assumed 
that  this  signal  is  digitized  and  converted  into  a  floating  point  number.  In  this  thesis 
it  is  assumed  that  the  floating  point  arithmetic  is  of  sufficient  precision  that  rounding 
errors  are  not  a  significant  problem. 

B.     PROBLEM  DEFINITION 

The  signal  described  in  section  (I.  A. 3)  can  be  thought  of  as  a  broad-band  signal 
with  additive  narrow-band  signals  occurring  at  arbitrary  frequencies  and  times.  These 
narrow-band  signals  can  be  modeled  as  either  pure  sinusoidal  signals,  as  outputs  of  a 
narrow-band  filter  excited  by  white  noise,  or  as  any  band-limited  narrow-band  signal. 

We  shall  assume  that  the  signal  y^  is  received  by  the  receiver  and  is  defined  as 

3 

yk  =  wk  +  ^2x\  +  ^k  (l.l) 

where  xj-.,  the  interferer,  is  either  a  pure  sinusoid  or  a  narrow-band  signal  whose 
frequency  (i.e  center  frequency)  can  vary  slowly  with  time.  The  term  Wk  is  the 
broad-band  DSSS  signal  of  interest  and  7^  is  White  Gaussian  Noise(WGN)  with  zero 
mean  and  variance  a2  Ar(0,  a2)  The  goal  is  to  remove  only  the  narrow-band  interferers 


and  obtain  the  output  (wk  +  7*)  for  further  processing.  A  typical  spectrum  for  y^  is 
given  in  Fig  1.2. 

In  this  thesis,  the  design  for  an  adaptive  digital  filter  that  will  achieve  the  goal 
of  removing  narrow-band  interferers  from  a  broad-band  signal  of  interest  is  developed. 
This  thesis  is  organized  into  five  chapters,  chapter  I  is  the  introduction.  The  chapter 
II  concentrates  on  building  necessary  background  on  such  things  as  the  Least  Square 
Estimator  and  Maximum-Likelihood  Estimator.  Chapter  III  briefly  explains  existing 
methods  of  adaptive  filtering  with  emphasis  on  adaptive  notch  filtering  techniques 
.  The  new  algorithm  [Ref.  26]  developed  for  the  purpose  of  addressing  the  current 
problem  of  interest  is  also  explained.  Chapter  IV  is  totally  devoted  to  modelling 
various  signals  required  for  testing.  Simulation  results  are  discussed  and  presented. 
Chapter  V  provides  a  brief  survey  of  existing  hardware  components  and  a  look  at 
possible  hardware  structures  and  their  limitations.  Chapter  VI  is  conclusions. 
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II.  ESTIMATION  TECHNIQUES 

A.     DIGITAL  FILTERS 

Digital  Filters  can  be  characterized  by  their  impulse  response,  their  trans- 
fer function,  or  by  a  difference  equation  [Ref.  13]  .  A  typical  Infinite-Impulse- 
Response(IIR)  filter  represented  in  difference  equation  form  is  given  by 

xk  =  aizjt-i  +  a2xk-2  +  a3Xfc_3  +  bxuk  +  62u*-i  (2.1) 

while  an  Finite  Impulse  Response(FIR)  filter  is  given  by 

Xk  =  biiik  -f  &2Ujt-i  +  b3uk-3  (2.2) 

where  uk  is  the  input  to  the  filter  and  xk  is  the  output  of  the  filter.  xk  and  uk  are 
the  sampled  values  of  the  continuous  functions  x(t)  and  u(t).  By  choosing  Ta  as  the 
sampling  time  we  notionally  write  the  discrete  signal  as 

xk  =  x{k)  =  x(kT3)  (2.3) 

Let  X(Sl)  be  the  spectrum  of  the  signal  xk,  where  fi  is  the  angular  frequency 
component  present  in  the  actual  signal.  For  analysis  purposes  consider  normalized 
frequency  given  asw  =  ClTa.  With  this  definition  we  do  not  have  to  use  the  actual  fre- 
quencies but  only  the  normalized  frequencies.  Assuming  that  there  is  no  aliasing  while 
sampling  the  maximum  normalized  frequency  is  0.5  cycle/sample  or  ir  rad/sample. 

1.      Some  Hardware  Schemes 

Fig  2.1  represents  a  conventional  Direct  Form  II  representation  of  a  digital 
filter  [Ref.  18].  However  for  high-speed  throughput  it  is  often  necessary  to  have  a 
pipelined  architecture.  FIR  filters  are  highly  amenable  for  pipelining  thus  giving  a 
high  throughput  rate  for  computations. 


2.  Modelling  Hardware  Multiplier  and  Adder 

Multiplication  of  two  variables  say  61  and  uk  in  hardware  can  be  modelled 
as 

bi*uh  -  hi  *  ukDx  =  bxukDi  (2.4) 

where  *  is  the  hardware  multiplier  and  *  is  the  ideal  multiplier  and  D\  is  the  propa- 
gation delay  associated  with  the  multiplication  Similarly  a  hardware  adder  is  desig- 
nated as  -f  while  the  ideal  adder  is  given  as  +  and  the  corresponding  delay  associated 
with  the  adder  is  D2.  Then  an  addition  of  two  variables  uk  and  uk-\  can  be  written 
as 

uk+uk-i  =  (uk  +  uk-\)D2  (2.5) 

Each  of  the  previous  devices  can  be  incorporated  into  a  pipeline  digital  filter  structure 
by  following  it  with  a  synchronizing  register.  The  clock  period  must  then  be  less  than 
D  4-  ts  +  tr  where  ts  is  the  register  setup  time  and  tr  the  register  propagation  delay. 

3.  Pipelining  an  FIR  filter 

In  this  section  we  wish  to  realize 

xk  =  6j  *  uk  +  62  *  tifc-i  +  63  *  Ujt_2  (2.6) 

It  is  acceptable  to  have  an  overall  delay  DN  ,  However  it  is  desirable  to  minimize  TV: 

xk  =  DN  {61  *  uk  +  62  *  tu_i  +  h  *  uk-2]  (2.7) 

First  associate  one  D  with  each  multiplier  and  form  real  multipliers  by  the  notation 
D*  =  *.  Then  we  get 

xk  =  DN~l  {bxD  *  uk  +  b2D  *uk  -  1  +  b3D  *  uk.2)  (2.8) 

Now  form  the  real  multiplier 

xk  =  DN~l  {b.iuk  +  b2*uk  -  1  +  b3*uk-2}  (2.9) 


Now  we  use  one  additional  D  to  associate  with  the  adding  of  62*iu-i  and  63*u^._2. 
(Note:   It  is  important  to  choose  the  most  delayed  values  first  in  order  to  minimize 

M 

xk  =  DN~2  {Dbx*uk  +  £(62*ufc-i  +  b3iuk.2)}  (2.10) 

Now  form  the  real  adder 

Xk  =  DN~2  [Dbxiuk  +  (&2*Ujfe-l4-&3*Ujfe-2)}  (2.11) 

Now  we  use  an  additional  delay  D  to  associate  with  the  final  adder  to  get  a  real 
adder: 

xk  =  DN~3  {j96i*wJfc+(62*tzfe_1+63*Hfc-2)}  (2.12) 

Finally  we  choose  the  system  delay  such  that  D  =  z~l  and  replace  i^_2  = 
z~2Uk  and  uk-\  =  z~1Uk,  thus  we  get 

xk  =  z"A'+3  {z-'b^+iz-'biiu.+z-'hiu,)}  (2.13) 

Now  factor  z~l  out  of  the  expression  to  obtain 

**  =  z-N+2{b1iuh^{b2*uk^z-1b3*uk)}  (2.14) 

Note  that  TV  =  2  is  the  minimum  possible  value  of  A^  for  a  causal  filter.    Choose 
TV  =  2  then 

xk  =  {6i*ufc+(62*ujk+-"163*W)fc)}  (2.15) 

Equation  (2.15)  is  very  convenient  for  pipelining  and  is  given  one  to  one  in  Fig  2.2. 
The  above  procedure  is  general  for  any  FIR  filter.  Pipelining  an  FIR  filter 
results  in  maximizing  the  sampling  or  throughput  frequency  at  the  expense  of  a  delay 
or  latency  in  the  availability  of  the  output. 
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4.      Pipelining  IIR  filters 

Pipelining  an  IIR  filter  is  considerably  more  difficult  than  an  FIR  filter. 
This  is  mainly  because  the  delay  in  the  availability  of  the  output  makes  it  impossible 
to  feed  back  output  values  with  short  delay  as  is  required  for  the  straightforward  IIR 
implementation.  As  an  example  consider  a  simple  1st  order  IIR  filter 

xk  =  ax  *  Xk-i  +  Uk  (2.16) 

As  in  the  FIR  case  we  introduce  a  delay  Dh  as  follows. 

xk^DN{a1*xk.1+uk)  (2.17) 

Now  we  use  one  delay  to  form  the  real  multiplier: 

xh    =    DN-1{alD*xk.1  +  Duk)  (2.18) 

=    DN-1{a1*xk.1  +  Duk)  (2.19) 

We  use  one  more  D  to  form  the  real  adder 

xk  =  DN-2{alD*xk.l+uk)  (2.20) 

Finally  assume  D  =  z~l  and  xk-\  =  z~lxk  then 

xk     =    z-^ia^z-'xk+z-'u,)  (2.21) 

=    z-N+1(aiixk+uk)  (2.22) 

For  xk  =  xk,  N  must  be  zero.  However  TV  >  1  is  required  for  a  causal  filter. 

Solutions  to  this  problem  exist.  They  use  a  higher  order  difference  equation 
representation  of  the  filter  with  equivalent  characteristics,  so  that  feedback  loop  delay 
greater  than  1  can  be  tolerated.  More  details  can  be  found  in  the  literature  [Ref.  14] 
regarding  pipelining  of  IIR  digital  filters. 
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B.     SIGNAL  CHARACTERIZATION 

The  signal  yk  of  equation  (1.1)  can  be  characterized  either  in  the  time  domain 
or  frequency  domain  and  the  signal  can  be  easily  transformed  from  one  domain  to  the 
other  via  a  linear  transformation. 

1.      Time  Domain  Representation 

In  the  time  domain  representation,  the  signal  yk  is  modelled  as  the  output 
of  a  linear  system  which  may  be  either  an  all  pole  or  a  pole-zero  system  [Ref.  13]. 
The  input  is  assumed  to  be  white  noise  but  this  input  to  the  system  is  not  accessible 
and  is  only  conceptual  in  nature.  Let  the  system  under  consideration  be 


V  9 


Vk  =  Yl  a*'V*-«  +  H  &i+i7fc-j  (2-23) 

t=l  j=0 

where  7^  is  a  white  noise  and  yk  is  the  output  of  the  system.    Equation  (2.23)  also 
can  be  represented  as  a  Transfer  Function (TF)  [Ref.  17] 

*M  =  fg  (2.24) 


or 


**  "l^rlk  (2-25) 


A(z) 

where 

p 
A(z)  =  1-  J^aiZ-*  (2.26) 

1=1 

and 


£(*)  =  £  W  (2.27) 

3=0 

we  define  the  parameter  vector  as 

p'  =  [ai,a>2,  •••»aj>»  hi  b2,-  —  ,bq]  (2.28) 

The  parameter  p  completely  characterizes  the  signal  Xk  The  above  system  defined 
by  the  equation  (2.23)  is  also  known  as  an  Auto-Regressive-Moving-Average(ARMA) 

12 


model  [Ref.  12],[Ref.  15].  It  is  conventional  to  denote  a  p  poles  and  q  zeros  system 
as  ARMA(p.q). 

2.      Frequency  Domain  representation 

The  signal  yk  of  equation  (2.23)  can  also  be  represented  as  a  vector 

YN  =  [yi..-yN]  (2.29) 

which  can  be  transformed  into  the  frequency  domain  by  the  DFT  relation 

»-^{S  ("""V-)}  (2-30) 

where  W  =  e-^  and  j  =  \f—\.  In  general  gk  is  a  complex  quantity.  However  we 
restrict  our  interest  to  the  power  spectrum  a  real  quantity  given  as 

Sk  =  9k9l  (2-31) 

The  series  Sat 

S,v  =  [*i  •••**]  (2.32) 

is  another  form  of  representing  the  signal  and  is  designated  as  the  power  spectrum  of 
the  signal  Xk-  Even  though  the  signal  in  this  domain  is  very  convenient  to  handle, 
computational  complexity  and  other  considerations  limit  its  use  for  online  application. 
Also  in  this  representation  phase  information  of  the  signal  is  lost. 

C.      SPECTRUM  ESTIMATION 
1.      Least  Square  Estimator- 
Consider  the  signal  yk  the  output  of  a  linear  system  excited  by  a  white 
noise  as  given  by  the  equation  (2.23).  Now  our  problem  is  to  estimate  the  parameter 
vector  p  by  obtaining  successive  measurements  of  yk-  We  shall  define  the  vector 

*l  =  btk-u-  •  •  *Vk-p,1[k,-  •  •  tfk-q]  (2.33) 

13 


and  Zk  =  Vk  so  that  the  difference  equation  (2.23)  can  be  rewritten  as 

zk  =  pUk  (2.34) 

an  elegant  solution  [Ref.  2]  for  the  above  problem  can  be  written  as 

Pk  =  Pfc-i  +  H^e*  (2.35) 

where 

h,  =  fe>txn  (2-36) 

and 


ejfc  =  zk-  Pfc-iXi  (2.37) 

The  block  solution  for  the  equation  (2.34)  can  be  also  written  in  the  form 

p*  =  Hfcfex*z*J  (2.38) 

Equation  (2.37)  can  also  viewed  as  Axk  —  B^k  and  this  is  represented  in  Fig  2.3. 

Computing  the  matrix  H^  via  equation  (2.36)  is  computationally  inefficient 
and  a  recursive  solution  via  the  matrix  inversion  lema  [Ref.  3]  is  preferred  and  is  given 
by 


Hfc  =  H^_i 


H/c_1x*xj.H/t_] 


(2.39) 


1  +x<fcH*_1xit  _ 
Staticians  refer  to  the  H^  matrix  as  the  covariance  matrix  while  optimization  people 

refer  to  H^  as  the  Hessian  of  the  objective  function,  where  the  objective  function  is 

given  by 

Jk  =  Y.e]  (2.40) 

Appendix  A  gives  a  computer  program  based  on  this  approach  that  was  used  for  this 
thesis  work  to  identify  system  parameters  via  Equations  (2.35),  (2.36),  (2.37),  (2.39). 
The  above  method  is  known  to  work  well  as  long  as  the  measured  value  of  the  state 
of  the  filter  is  free  from  noise.  In  the  event  zk  is  not  noise  free,  estimates  are  known 
to  be  biased. 
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2.      Maximum  Likelihood(ML)  Estimator 

The  ML  estimator  for  the  above  problem  is  conventionally  obtained  by 
considering  a  log-likelihood  function  [Ref.  16]  and  minimizing  it.  However  for  getting 
a  better  physical  understanding  the  approach  given  by  Young[Ref.  2]  will  be  followed. 
Consider  a  system  as 


Vk  = 


B(z) 


Ik 


and  an  auxiliary  system  referred  to  as  the  model  is  given  as 

B{z) 


(2.41) 


xk  = 


A(z)l 


Ik 


(2.42) 


and  we  define  the  function  to  be  minimized  as  L  =  e£(p)  where  ek  =  yk  —  xk. 
To  minimize  this  function  the  derivatives  and  the  Hessian  of  the  function  L  where 
L  =  (yk  —  xk)  need  to  be  obtained.  First  take  partial  derivative  of  L  with  respect  to 
each  parameter.  For  example 


—  =2(yk-xk)(-—)  =  -2ek(  —  ) 
The  partial  -^  can  be  obtained  by  differentiating  equation(2.42)  as 


(2.43) 


d    [B{z)     \  _d_ 

d&iU(z)7*J  dbx 


bi  +  b2z 


-l 


1 


U(*) 


1  —  a,\z~x  —  ci2Z~2  —  a^z~3 

Ik 


Ik 


Similarly  the  other  partial  -^  can  be  obtained  as 

8    \B{z)     \  ( A{z)^{B(z)}         B(z)z-* 

b(z)  _; 

■z      l  7^ 


\A(z)]' 


(2.44) 
(2.45) 


(2.46) 

(2.47) 
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since  -J^-{B(z)}  =  0  as  B(z)  has  no  terms  containing  cl\.  Using  equations  (2.45)  and 
(2.47)  the  partials  are: 


d 


Xk 


dik 
dyk 


=       Ul.   = 


=      II  = 


=    Vk  = 


1 


B(z) 


Ik 


i[A(zW 
1 


7fc 


M*)l 


Vk 


by  defining 


Pi     =     [ai»a2, 03,61,62) 

Xk      =      lxl-lixl-2ixl-3i  uliul-l] 

4  =  [y*-i,y)fe-2>y]fe-3>u*>wA-i] 

the  Hessian  may  be  approximated  by: 

Using  the  above  equations  the  parameter  can  be  estimated  via 

pk  =  Pk-i  +  ^Hjtx*e;t 


(2.48) 
(2.49) 

(2.50) 

(2.51) 


(2.52) 
(2.53) 
(2.54) 


(2.55) 


(2.56) 


where  fi  is  the  step  size  in  incrementing  the  parameter  vector.  A  block  schematic  for 
the  algorithm  was  given  in  Fig  2.4.  A  serious  drawback  of  this  method  is  the  stability 
while  incrementing  the  parameters. 
3.      Other  Methods 

There  are  many  other  methods  such  as  the  Instrumental  Variable  (IV) 
approach,  the  Approximate-Maximum-Likelihood  (AML)  Method  and  a  combination 
method  IV- AML  available  in  the  literature  [Ref.  2]. 
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D.      SELECTION  OF  ESTIMATION  TECHNIQUE 

Before  we  select  the  appropriate  method  for  the  problem  of  interest,  model-order 
problems  must  be  considered. 

1.      Model  Order  Problems 

In  the  Least-Square  Estimator  for  the  parameter  vector  p  of  the  signal 
Xk  we  need  to  assume  the  dimension  of  p  or  the  order  of  the  system.  Since  the 
system  order  is  often  not  known  a  higher-order  model  estimate  is  assumed.  Now  we 
investigate  higher-order  model  fits  for  a  lower-order  data. 

Consider  the  system  equation  (2.23)  with  parameter  given  as 

p*  =  [1.7,-1.53,0.648,1,0.6]  (2.57) 

For  the  above  system  a  Random  Binary  Noise(RBN)  was  given  as  the  input  Uk  and 
the  corresponding  output  Xk  is  obtained.  For  the  time  series  it*  and  Xk  an  ARMA(6,1) 
model  was  fit  using  the  program  given  in  the  appendix  and  the  parameter  was  esti- 
mated as 

p'  =  [2.26533,-2.79577,2.17249,-1.09877,0.446598,-0.110151,0.998175]     (2.58) 

Fig  2.5  gives  the  parametric  spectrum  for  the  above. 

For  the  same  data  an  ARMA(4,1)  model  was  used  and  generating  a  pa- 
rameter vector 

p'  =  [2.16988,-2.41286,1.47436,-0.365412,0.999317]  (2.59) 

Fig  2.6  gives  the  parametric  spectrum  for  the  above  values. 

Lastly  for  the  same  data  an  ARMA(8,1)  model  was  used  and  generating  a 
parameter  vector 

p'  =  [2.276,  -2.835, 2.267,-1.282, 0.703,  -0.353, 0.143,  -0.0035, 0.9999]        (2.60) 
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Fig  2.7  gives  the  corresponding  spectrum. 

It  is  interesting  to  compare  these  with  the  actual  parametric  spectrum  given 
in  Fig  2.8  For  the  sake  of  completion  and  comparison,  the  spectrum  given  in  Fig  2.9 
of  the  output  time  series  Xk  computed  via  a  DFT  program  is  also  included.  This 
simulation  demonstrates  that  the  existence  of  multiple  solutions,  hence  an  important 
step  in  the  estimation  procedure  is  to  choose  an  appropriate  model  order. 
2.      Choice  of  the  Method 

Since  the  filter  hardware  must  be  implemented  in  real-time  at  the  frequen- 
cies of  10  to  20  Mhz  computational  complexities  must  be  kept  to  a  minimum.  These 
methods  although  providing  good  performance,  are  not  well  suited  to  hardware  im- 
plementations. This  motivates  looking  into  new  filtering  methods. 
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III.  METHODS  for  FILTERING 

A.     FILTERING  TECHNIQUES 

Filtering  in  a  broad  sense  is  selectively  suppressing  a  portion  of  the  spectrum 
of  the  given  signal.  In  this  chapter  several  filtering  techniques  are  explored  in  the 
light  of  the  narrow-band  interference  problem  in  order  to  identify  applicable  filtering 
techniques. 

1.      Filtering  via  FFT 

Any  given  signal  can  be  conveniently  transformed  into  its  power  spectrum 
via  equation  2.31.  Assuming  the  spectrum  of  the  desired  filter  to  be  defined  by  w^. 
The  weighted  output  is  given  by: 

vk  =  skwk  (3.1) 

and  the  corresponding  filtered  output  y*  is  obtained  by  taking  the  inverse  DFT  of  the 
signal  Vk-  A  simple  block  schematic  representing  this  idea  is  given  in  Fig  3.1.  Details 
of  this  approach  are  available  in  various  references  [Ref.  22]. 

Fast  algorithms  such  as  the  FFT  for  computing  the  DFT  make  it  possible 
to  do  the  above  process  in  real  time  by  using  dedicated  hardware.  Honeywell  makes 
the  HDSP66110  and  HDSP66210  Digital  Signal  Processing  chip  pair  which  are  ideally 
suited  for  these  applications.  This  chip  pair  can  perform  a  single  complex  multiply  in 
40ns  [Ref.  28].  However  the  problem  of  filtering  adaptively  [Ref.  10]  (i.e  varying  Wk 
according  to  a  criterion)  demands  more  computing  power  which  can  be  obtained  by 
augmenting  the  DSP  chip  pair  with  a  processor  like  the  R3000  which  has  a  Reduced 
Instruction  Set  Computer(RISC)  architecture  with  a  high  instruction  execution  rate 
of  approximately  25  Million  Instructions  Per  Second(MIPS)  [Ref.  27]. 

23 


2.  Recursive  DFT 

Implementation  of  a  higher  order  FIR  filtei  using  a  recursive  I > I  1   [Ref 
23], [Kef.   24]  is  also  very  convenient  foi  eliminating  narrowband  interference.    The 
basic  concept  is  thai  the  FIR  filtei  is  expressed  as  a  producl  of  two  filtei  sections 

One  section  is  a  filter  with  its  zeros  being  equally  spaced  on  the  unit  circle.  I  his 
is  achieved  by  a  delay.  The  second  section  is  a  pole-producing  section.  Pole  zero 
cancellation  results  in  the  desired  FIK  filter.  More  details  can  be  seen  in  the  referen<  es 
[Ref.  23], [Ref.  24]. 

3.  Adaptive  Filtering 

Adaptive  filters  can  be  placed  into  fom  classes  based  upon  the  choice  of 
the  training  sequence  and  the  reference  model  fol  adaptation.  Simon  \\{<-[.  21]has 
classified  the  Adaptive  Filters  into  the  following  four  classes 

•  Identifical ion  -  ('lass  I 

•  Inverse  Modelling  -  Class  II 

•  Prediction  -  Class  III 

•  Interference  Canceling  -  Class  IV 

Fig  [3.2  -  3.5]  give  the  block  schematics  of  the  various  classes  of  adaptive 
filters.  In  oui  problem,  the  reference  model  is  naturally  a  narrow-band  bandpass  filtei 
since  the  interference  signal  is  a  narrow -band  signal.  The  class  of  filtering  thai  besl 
suits  our  problem  is  a  combination  of  Class  III  and  Class  IV  type  of  filtei  [Fig  3.4 
and  3.5].  I  his  logically  points  to  an  adaptive  notch  filter.  Due  to  apnori  knowledge 
of  the  interference  signal,  only  the  parametric  approach  was.  However  the  DFT- 
based  techniques  may  also  offer  a  good  solution  and  should  be  the  subje<  I  of  future 
investigal  ions. 
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B.     ADAPTIVE  NOTCH  FILTERS 

Notch  filters  for  removing  multiple  narrow-band  interference  can  be  categorized 
into  four  broad  categories  illustrated  in  Fig  3.6  through  3.9.  The  first  two  categories, 
Figures  3.6  and  3.7,  are  cascaded  second  order  notches  with  each  second-order  section 
removing  one  frequency.  The  next  two  categories, 

Figures  3.8  and  3.9,  are  higher-order  notches  that  eliminate  multiple  frequen- 
cies. In  all  of  the  categories,  it  is  possible  to  use  FIR  filters  (i.e  all  zero  filters)  which 
are  easily  pipelined  and  can  be  made  truly  linear  phase.  However,  IIR  notch  filters 
require  substantially  fewer  multipliers  and  adders  than  FIR  notch  filters.  Thus  IIR 
pipelining  may  become  an  important  issue.  In  this  thesis  we  limit  our  discussion  only 
to  the  first  two  categories  illustrated  in  Figures  3.6  and  3.7. 

1.      Second-Order  Cascaded  Notch  Filters 

The  second-order  notch  filter  is  used  in  cascade  and  in-line  with  the  signal 
as  shown  in  Figure  3.6.  The  transfer  function  for  such  a  notch  filter  is  given  by  Kwan 
and  Martin  [Ref.  8]  as: 

HN(z)    =    l-Hbp(z)  (3.2) 

«    l-I*i (1 +  «-»)(!-*-') )    ■ 

\  2  1  -  (2  -  k2  -  k\)z-i  +  (1  -  k2)z-*  J  K      } 


2-h  1  -  ^j^-1  +  z~2 


(3.4) 


2  1  -  (2  -  k2  -  k\)z-1  +  (1  -  k2)z~2 
For  arbitrary  values  of  k^  and  k2  ,  this  is  a  symmetric  notch  filter  with  unity 
gain  at  DC  and  the  Nyquist  frequency.  If  k2  is  kept  constant,  then  the  3db  notch 
width  is  also  kept  constant.  Thus  ki  may  be  adapted  to  remove  one  narrow-band 
signal.  A  cascade  of  such  filters  can  be  used  to  remove  multiple  narrow-  band  signals 
[Ref.  8].  Constants  kj  and  k2  are  related  to  the  pole  radius  r  and  the  normalized  pole 
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frequency  6P  as 

k2    =    (1-r2)  (3.5) 


k,     =     y/\  +  r2  -  2rcos6p  (3.6) 

It  is  important  to  bear  in  mind  that  the  above  transfer  function  has  unity  gain  at 

Opeak  =  2arcsin  <  — .   1         >  (3.7) 

l2v/n^j 

which  is  different  from  6P  [Ref.  8] 

2.      Second-Order  Cascaded  Signal  Canceler 

The  cascaded  second-order  signal  canceler  approach  shown  in  Figure  3.7 
has  the  advantage  that  the  desired  signal  does  not  pass  through  the  adaptive  filter. 
Instead,  the  band-pass  filter  is  used  to  detect  the  narrow-band  signal  which  is  then 
subtracted  from  the  desired  signal.  A  constant  3db  bandwidth  notch  can  be  achieved 
by  selecting  a  band-pass  filter  with  the  transfer  function: 

*M    =     "-^P  (3-8) 

-     £&  (39) 

D(z)  (39) 

where  D(z)  is  the  same  denominator  as  H^(z)  in  equation  (3.4)  The  signal-canceler 
structure  is  also  nice  for  adaptation  because  it  is  relatively  easy  to  generate  sensitivity 
functions  which  are  related  to  the  gradient  of  Hbp{z)  with  respect  to  the  frequency 
parameter  k\.  Figure  3.10  shows  the  block  diagram  of  an  adaptive  version  of  this 
filter.  The  sensitivity  function  s(n)  is  obtained  by  differentiating  the  error  signal 
e(n)  with  respect  to  the  parameter  k\.  This  can  be  easily  derived  by  noting  that 
e(n)  =  x(n)  —  y(n)  and 

oki  oki 
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to  get  the  above  partial  derivative  we  consider 

y(n)    =    Hbp(z)x(n) 
N(z) 


D(z) 


x(n) 


(3.11) 
(3.12) 


the  above  equation  (3.12)  can  be  easily  differentiated  with  respect  to  ki  by  using  the 
equation  (2.47)  to  obtain: 


dy(n)  =       [  N(z) 


2kiz 


-i 


x(n) 


(3.13) 


by  recognizing  that  Hbp{z)  =  ^jfj  from  equation  (3.12)  the  above  equation  can  be 
written  as 


dyjn) 


Hbpiz) 


2kxz 


-i 


and  by  defining 


Has(z)  = 


D(z) 
2k1z~1 


x(n) 


(3.14) 


(3.15) 


we  get  the  sensitivity  function  as 


s(n)  =  ?jj£l  =  H.(z)  =  Hbp(z)HS3(z) 


(3.16) 


The  equation  (3.16)  is  given  as  a  block  schematic  in  Fig  3.10.  The  parameter  k^  may 
then  be  adapted  by  the  formula: 


S\  Tl ) 

kx{n  +  1)  =  h(n)  -  ne{n)  2 

ll5(")ll 

a.       Computing  ||s(n)||2 

Ideally  ||5(n)||    can  be  calculated  by  the  relation 


(3.17) 


\s{n)\\2=    £    s(i)2 

i=n-N 


(3.18) 
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The  above  scheme  computes  the  average  of  the  past  N  samples.  However  a  weighted 
average  of  the  past  samples  with  highest  weight  for  the  recent  sample  can  be  done 
quite  easily  by  a  first-order  low  pass  filtering  [s(n)]2  as  follows: 

vn  =  vn-1\  +  {l-  X)s2{n)  (3.19) 

However  the  above  lowpass  filtering  can  also  be  carried  out  by  a  second  order  filter 
such  as: 

„„  =  *»('.- 2r'.+z7  mr  (3.20) 


1  -  2rfkz~l  +  z 


vn  = 


-2 
-1   j     ~-2- 


6i(l  -  2cos{59)z-1  +2-2) 


1  -2rjkz~l  +  z~2 
Yet  another  way  to  estimate  ||s(n)||    i  simply: 

N 


[s(n)]2  (3.21) 


»»  =  2>W  (3-22) 


the  above  form  is  defined  as  zero  order  forgetting.  It  should  be  noted  that  equation 
(3.21)  was  used  by  Kwan  and  Martin  [Ref.  8]. 

C.      MULTIPLE  NARROWBAND  SIGNAL  SUPPRESSION 

The  second-order  implementations  of  section  (1.)  and  (2.)  offer  considerable 
advantage  both  in  hardware  complexity  and  in  adaptive  performance  for  multiple 
narrow  band  signal  suppression  [Ref.  8], [Ref.  25], [Ref.  26]. 

1.      The  Kwan  and  Martin  Filter 

In  a  recent  paper  by  Kwan  and  Martin  [Ref.  8],  the  problem  of  detecting 
and  enhancing  sinusoidal  signals  in  the  presence  of  noise  is  addressed  with  a  cascade 
of  IIR  adaptive  notch  filters  which  are  used  to  eliminate  the  sinusoids.  Each  of  the 
sinusoids  is  eliminated  by  a  bandpass  filter  whose  output  is  an  enhanced  version  of 
one  of  the  sinusoids.  Hence  this  remarkable  structure  can  perform  both  tasks  with  a 
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single  adaptive  filter  configuration  which  is  shown  to  be  highly  robust  and  performs 
extremely  well. 

The  major  disadvantage  of  the  Kwan  and  Martin  structure  is  that  the 
number  of  biquad  sections  needed  in  the  adaptive  filter  configuration  is  given  by 
N(N+3)/2,  where  N  is  the  number  of  sinusoids  to  be  detected  and  removed.  This 
becomes  impractical  in  real-time  situations  with  more  than  4  sinusoids  due  to  the 
geometric  increase  in  the  required  hardware. 

a.      Kwan  and  Martin  Structure 

The  Kwan  and  Martin  structure  consists  of  a  cascade  of  IIR  notch 
filters  one  stage  of  which  is  shown  in  Figure  3.11.  Each  stage  consists  of  a  bandpass 
filter  with  zeros  at  DC  and  the  Nyquist  frequency  and  unity  gain  at  its  peak  frequency 
u>i.  Such  a  filter  would  have  the  following  z-domain  transfer  function 


HlP  =  -^-^ o a  —i    ,    ..2—7  (3-23) 


1  -  r}  1  -  z-' 

2      1  -  2rtcos6xz-     \-  r^z~2 

where 


r,    =  pole  radius  of  the  i-th  section 
6t  =    2"KUi/us 

lu\  =    peak  frequency  of  the  i-th  section 
u;s  —    sampling  frequency 

Kwan  and  Martin  identify  two  different  methods  for  adapting  the 
filter.  Most  of  their  derivation  is  based  on  what  they  call  the  constant  bandwidth  filter 
in  which  the  pole  radius  r,  is  a  constant  and  only  the  frequency  u>,  is  adapted.  An 
alternative  approach  which  keeps  a  constant  Q  is  also  discussed  in  Kwan  and  Martin 
paper.  In  addition,  Kwan  and  Martin  select  the  adaptive  quantity  in  such  a  way  that 
it  is  fairly  easy  to  determine  the  notch  frequency  from  the  adaptive  parameter.  From 
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Figure  3.10,  we  see  that  the  notch  filter  for  each  section  is  the  difference  between  1 
and  the  bandpass  filter,  hence: 

H*N{z)  =  1  -  Hlp(z)  (3.24) 

b.      Calculation  of  the  gradient 

The  basic  structure  of  the  Kwan  and  Martin  adaptive  filter  shown  in 
Figure  3.12  is  a  cascade  of  N  sections  of  the  form  of  Figure  3.7.  The  overall  transfer 
function  is  given  by 

T{z)    =    UH>N(z)  (3.25) 

t=i 

=    fl(l-0  (3-26) 

t=l 

Kwan  and  Martin  choose  as  their  objective  function  J(z)  the  square 

of  the  output  of  the  final  stage  of  the  cascade: 


J(z)    =    {E1(z)}2  (3.27) 

=    {T(z)]2[X(z)}2  (3.28) 

Hence,  the  gradient  of  the  objective  function  J(z)  is  given  by 

3gl -«(.„««£!  (3,9) 

Thus  in  order  to  find  the  gradient  with  respect  to  the  adaptive  pa- 
rameters kj  ,  we  must  take  the  partial  derivative  of  T(z)  with  respect  to  each  kj 

-qt-  -  Il^CO-gr-  (3-3°) 
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From  equation  (3.24)  we  have 

dHfo)     =     g[l  -  Hjp(z)} 
dkj  dkj 

dHjv(z) 


(3.31) 
(3.32) 


Using  the  rule  for  differentiation  as  given  by  the  equations  (2.45)  and  (2.47)  and  from 
equation  (3.23)  we  have 


dH},(z)  •       -Ikjz 

=    H3hp{z)HUz) 


-i 


dkj 


(3.33) 
(3.34) 


where  D(z)  is  the  denominator  of  HJhp.  HJa3  is  given  by  the  equation  (3.15). 

Substituting  equation  (3.34)  into  equation  (3.30)  we  obtain  the  over 
all  sensitivity  for  the  jth  ki  parameter  as 


=    UH'y(z)Hlp(z)Hi3(z) 


dk, 


;'-2 

N 

n^-(^) 

n  m*) 

t=i 

_»=J-1 

Htiz)Hi,{z) 


By  recognizing  that 


n  mz)Htt 


br 


Y>(z)  = 

and  by  substituting  equation  (3.37)  in  the  equation  (3.36)  we  get 

dT(z) 


(3.35) 
(3.36) 

(3.37) 


dk. 


H>3(z) 


(3.38) 


Figure  3.12  shows  the  Kwan  and  Martin  realization  of  the  complete 
adaptive  system.  The  difficulty  in  the  Kwan  and  Martin  approach  is  in  the  generating 
of  the  product  of  notch  filters  without  the  notch  filter  j,  as  required  in  equation  (3.36). 
To  generate  this  product  for  each  section,  would  require  TV  —  1  biquads  per  section 
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resulting  in  a  total  of  Ar(Ar  —  1)  biquads  just  to  generate  the  product.  Kwan  and 
Martin  are  able  to  reduce  this  by  using  the  output  of  the  bandpass  filter  as  the  input 
to  their  cascade  via  equation  (3.37).  Since  this  output  already  has  (j  —  1)  of  the 
required  Hfr(z)  factors  in  it,  only  (A1  —  j)  additional  biquads  are  needed  for  a  total  of 
0.5(A'2  —  A7).  Adding  this  to  the  N  biquads  required  to  realize  the  cascade  of  notch 
filters  and  the  Ar  biquads  required  to  realize  the  H^(z)  factors,  yields  a  total  number 
of  biquads  given  as  0.5(Ar  -j-  3)Ar. 
2.      The  New  Structure 

Figure  3.13  shows  the  improved  adaptive  notch  filter  structure  [Ref.  25],[Ref. 
26].  The  key  to  the  improvement  is  the  recognition  that  the  output  E\(z)  =  T(z)X(z) 
for  the  cascade  of  the  notch  filters  can  be  written  both  as  a  product  of  the  individual 
notch  filter  section  transfer  functions  H*N(z)  times  the  input  and  in  terms  of  the  input 
X(z)  minus  the  outputs  Yl{z)  of  the  bandpass  filters 


E,{z)  T(z)X(z) 


t=i 


X(z) 


=  A'(z)-(]TrW) 


(3.39) 
(3.40) 

(3.41) 


To  get  the  product  of  HxN(z)  without  the  term  i  =  j  as  required  in  the 
equation  (3.36),  we  may  use  equation  (3.41)  to  simply  add  back  the  term  YJ(z) 


iW(~~) 

X(z) 

=    X(z)- 

N 

1=1 

=    £1(2) + 

t=i 

(3.42) 
(3.43) 


Figure  3.13  makes  use  of  this  fact  to  generate  the  gradient  needed  for  the 
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adaptive  process.  From  the  Figure  3.13,  we  can  see  that  the  total  number  of  biquads 
required  is  N  for  the  cascade  of  notch  filters  plus  2N  for  the  Hxbp(z)H\{z)  required  for 
adaptation,  minus  1  at  the  last  stage,  since  the  last  stage  does  not  need  the  extra 
Hxb  (z)).  Thus  we  have  the  number  of  biquads  required  in  the  new  structure  is  (3N-1). 
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IV.  MODELLING  and  SIMULATION 

A.     MODELLING  OF  SIGNALS 

In  order  to  test  various  algorithms  and  to  evaluate  their  probable  performance 
in  the  real  environment,  it  was  necessary  to  develop  a  meaningful  simulation  of  the 
real  situation.  To  achieve  this,  the  following  four  testing  categories  were  developed: 

(i)  Sinusoidal  signals  with  white  gaussian  noise 
(ii)  Narrow  Band  Noise  with  white  gaussian  noise 
(iii)  Bi-Phase  Shift  Keying  (BPSK)  sequence 
(iv)  Frequency  Shift  Keying  (FSK)  sequence 

1.  Sinusoidal  Signals 

In  order  to  generate  sinusoidal  signals  with  minimum  computational  burden 
placed  at  different  frequencies  ^§5/*,  a  second  order  AR  process 

x\=2cos{9x)x\_,-xik_2  (4.1) 

was  used.  Initial  conditions  are  very  important  and  they  are  chosen  such  that  X-\  =  0 
and  X-2  =  —sin(Oi)  giving  a  unit  amplitude  sinusoidal  signal.  The  0,  value  is  between 
0  to  180  and  n  is  the  number  of  frequencies  desired.  The  required  signal  y^  needed 
to  input  into  the  adaptive  algorithm  is  given  as 

n 

lfr  =  X>i  +  1*  (4-2) 

t=l 

where  7^  is  a  white  gaussian  noise  7V(0,<72). 

2.  Narrow  Band  Noise 

The  narrow  band  noise  signal  is  generated  using  the  difference  equation 

z\  =  2rcos(0,)4-i  -  r2x\_2  +  uj  (4.3) 
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where  9t  decides  the  placement  of  the  noise  in  the  spectrum  and  r  controls  the  band- 
width of  the  noise.  The  u\  is  simply  a  uniformly  distributed  noise  taken  at  different 
instants.  The  desired  signal  yk  is  obtained  via 

Vks=I>i  +  7*  (4-4) 

t=i 

3.      Bi-Phase  Shift  Keying  Sequence 

The  generation  of  BPSK  signal  has  three  distinct  three  parts: 

(i)  Generation  of  Random  Binary  Sequence(RBS)  is  achieved  by  passing 

uniformly  distributed  noise  through  a  hardlimiter  (An  important  note  is  that  the 

interval  between  the  two  consecutive  bits  of  RBS  is  4-); 

lb' 

(ii)  Another  sequence  of  binary  numbers  is  a  spreading  code  or  sequence. 
The  specific  sequence  used  in  a  given  communications  system  is  normally  not  available 
to  anyone  but  to  the  designated  receiver.  (In  this  particular  simulation  we  have 
generated  the  spreading  sequence  by  passing  uniformly  distributed  noise  through  a 
hardlimiter.  The  bit  interval  is  4-); 

(iii)  Phase  encoding  (i.e.)  mapping  the  given  binary  signal  which  is  the 
exclusive  or  of  (i)  and  (ii)  as  0  or  tv  at  appropriate  sample  time. 

The  output  of  the  first  hardlimiter  is  stored  in  an  array  x.  Output  of  the 
second  hardlimiter  is  stored  in  the  array  y  This  information  is  retrieved  by  a  subtle 
use  of  the  array  index  given  as  i  =  kfb  where  i  is  the  index,  k  is  the  discrete  sample 
number  and  /{,  is  the  bit  rate  of  the  intelligence.  Similarly  another  index  j  is  generated 
using  j  =  kfc  where  fc  is  the  chip  frequency.  Generation  of  the  indices  is  the  key 
thing  in  this  simulation.  The  desired  signal  now  is  given  by  the  equations 

rvpk  =  Acos{2irfok  +  <j>k)  (4.5) 

4>k  =  MOffiyO)]*  (4.6) 
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p=n 


"*  =  IX  (4-7) 

P=i 

where  i  and  j  are  indices  of  the  arrays  as  defined  earlier.  It  is  an  important  point 
to  note  that  the  characteristic  of  the  signal  wvk  are  defined  by  the  parameter  p  = 

{/ch/cM- 

This  signal  is  not  really  a  simple  BPSK  signal  but  it  has  an  additional  fea- 
ture of  spreading  the  spectrum  by  controlling  the  chip-frequency  and  carrier  frequency 
and  information  rate.  A  block  diagram  of  the  scheme  is  given  in  Fig  4.1. 

The  desired  signal  yk  is  given  by: 

yk  =  ak  +  pk  (4.8) 

where  ak  is  the  set  of  narrowband  BPSK  signals  placed  at  different  places  in  the 
frequency  spectrum  and  (3k  is  the  broad  band  BPSK  signal  generated  for  a  specific  p 
value. 

4.      Frequency  Shift  Keying  Sequence 

Generating  this  sequence  needs  a  random  binary  intelligence  signal.  This 
was  once  again  is  achieved  by  passing  a  uniformly  distributed  noise  through  a  hardlim- 
iter.  The  output  of  this  hard  limiter  stored  in  an  array  x.  An  index  i  is  chosen  such 
that  i  =  kfb  where  /t  is  the  baud-rate  of  the  information  and  k  is  discrete  sample 
number.  Now  the  desired  signal  is  generated  via 

Sk  =  2cos(6k)sk-i  -  sk-2  (4.9) 

$k  =  0  +  6x(i)  (4.10) 

yh  =  Sk  +  ik  (4-11) 

where  6  is  the  carrier  frequency  and  6  is  the  depth  of  the  frequency  modulation. 
Initial  conditions  are  very  important  and  they  are  chosen  such  that  s_i  =  0  and 
5_2  =  —sin(B)  giving  an  unit  amplitude  sinusoidal  signal. 
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B.     SIMULATION 

The  adaptive  digital  filter  algorithm  as  described  in  the  earlier  chapters  is  sim- 
ulated and  tested  using  synthetic  data.  Simulation  was  carried  out  on  a  VAX  11/785 
computer  using  Fortran  77.  Listing  of  the  program  used  is  enclosed  in  the  appendix. 
The  adaptive  filter  parameters  of  interest  are 

(a)  Sharpness  of  the  notch  filter  defined  by  pole  position  (r2  =  1  —  n2) 

(b)  Step  size  in  the  parameter  update  procedure  (//) 

(c)  Time  constant  of  the  fading  filter  (A)  (  refer  to  equation  3.19) 

(d)  Model  order  {N=  #  of  2nd  order  filters) 

(e)  Order  of  the  incoming  signal  or  number  of  interferers  (m) 

1.      Kwan  &:  Martin  Algorithm 

In  simulating  this  algorithm,  the  most  important  thing  is  the  implemen- 
tation of  the  notch  filter.  Fig  4.2  gives  the  structure  of  the  algorithm.  Let  x(n)  be 
the  input  to  the  Adaptive  Filter  and  e(n)  the  desired  output.  This  desired  output  is 
obtained  by  passing  the  input  x(n)  via  a  cascade  of  TV  notch  filters  as 

e(n)  =  {H%(z)H»-*  •  •  •  Hh(z)}  x(n)  '  (4.12) 

Output  at  the  intermediate  jth  section  of  the  band-pass  filter  is  designated  as  y*(n) 
as  shown  in  Fig  4.2  Then  the  sensitivity  of  the  k{  parameter,  is  obtained  by  passing 
this  y-'(n)  through  another  cascade  of  (j  —  1)  notch  filters  given  as 

s>(n)  =  {H^Hj,-2  •  •  •  HhH'.}  y>(n)  (4.13) 

The  filter  transfer  functions  HJN(z)  and  Hjs  were  defined  in  the  earlier  chapter  by 
equations  (2.21)  and  (2.15). 
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2.  The  New  Algorithm 

The  primary  difference  between  the  New  Algorithm  and  the  Kwan  and 
Martin  algorithm  is  in  the  method  of  obtaining  the  sensitivity  of  the  jth  coefficient. 
This  difference  can  be  seen  in  the  Fig  4.3.  The  output  r/J(n)  at  the  jth  notch  Fig  4.3. 
filter  is  added  to  the  over-all  output  e(n)  and  this  summed  output  is  passed  through 
a  cascade  of  only  two  filters  as  shown  in  Fig  4.3  and  can  be  given  as 

s>(n)=  [H^HQ  {e(n)  +  yj(n)}  (4.14) 

3.  Forgetting  Filters 

For  parameter  incrementation  via  equation  (3.17)  we  need  to  obtain  ||s2(n)||. 
This  could  be  achieved  by 

a)  Zero  order  forgetting  using  equation  (3.22) 

b)  1st  order  forgetting  using  equation  (3.19) 

c)  2nd  order  forgetting  using  equation  (3.21) 

4.  Stability 

The  parameter  incrementation  given  by  the  equation  (3.17)  can  be  recast 
by  using  either  of  the  forgetting  filters  given  above  as 

K\  =  k\  -  //[t',]_1s,(n)e(n)  (4.15) 

It  is  very  important  to  note  that  the  above  equation  (4.15)  has  a  close  resemblance 
with  the  ML  Estimator.  In  the  MLE  case  the  stability  while  incrementing  the  param- 
eter is  an  important  factor.  A  similar  problem  exists  in  the  current  incrementation 
procedure  but  the  problem  is  solved  by  a  suitable  choice  of  parameters  and  modify- 
ing the  incrementation  procedure  by  adding  an  additional  factor  to  equation  (4.15) 
as  follows: 


k\  =  k\  —  fie(n)sx(n) 
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This  modification  also  protects  against  possible  overflow  or  underflow  while  computing 
[i>*]~  .  In  addition  simulations  indicate  that  in  the  case  of  zero-order  forgetting,  this 
factor  is  required  in  order  to  obtain  convergence.  When  i>m,-n  and  vmcLX  are  zero  and 
infinity  respectively,  equation  (4.16)  reduces  to  equation  (4.15). 

In  this  incrementation  procedure  two  forms  of  the  denominator  polynomial 
of  the  Ht,p(z)  were  considered: 

D\z)  =  1  -  (2  -  k2  -  k\)z-1  +  (1  -  k2)z~2  (4.17) 

and 

Dx{z)  =  1  -  2/cir*"1  +  r2z~2  (4.18) 

where  k\  =  cos(6'). 

While  using  the  incrementation  procedure,  under  transient  conditions  poles 
of  the  Dl(z)  cross  the  unit  circle  causing  instability  problems.  This  condition  was 
averted  by  checking  the  pole  position  after  incrementation  and  if  unstable  then  in- 
crementation is  modified.  Checking  this  condition  for  Dl(z)  given  by  equation  (4.17) 
calls  for  solving  a  quadratic  equation  at  every  parameter  increment  and  examining 
the  pole  position.  However  this  check  is  very  simple  for  Dx (z)  given  in  equation 
(4.18),  since  we  need  only  to  check  k\  by  maintaining  \k\\  <  1.  Furthermore  the 
value  of  r  must  be  positive  and  less  than  unity  for  stability.  The  choice  of  the  r  is 
very  important  for  proper  fast  convergence.  Fig  4.4.  shows  the  frequency  response  of 
band-pass  filter  for  different  values  of  r  The  3dB  band-width  of  the  band  pass  filter 
is  related  to  r  and  is  given  [Ref.  29]  by 


\i+W 


Band  -  width  =  cos'1  < ->  (4.19) 

Note  that  under  limiting  conditions,  the  band-width  is  essentially  zero.  In  this  simu- 
lation it  was  assumed  a  value  of  r  as  0.95. 
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A  large  number  of  simulations  [Ref.  29]  were  carried  out  using  sinusoidal 
inputs  with  and  without  noise  for  different  m  and  N  values.  In  these  simulations 
the  denominator  polynomial  used  was  given  by  equation  (4.18).  The  results  were 
tabulated  in  Table  3.1.  From  Table  3.1  it  is  seen  that  the  choice  of  the  step  size  is 
an  important  factor  and  step-size  must  be  optimized  for  a  specific  number  of  sections 
N.  The  values  of  t>mtn  and  vmax  did  not  pose  any  problem  while  incrementing  using 
1st  order  or  2nd  order  forgetting  filters.  It  seems  that  a  1st  order  forgetting  filter  is 
more  effective  than  either  zero-order  or  a  2nd  order  filter. 

After  fixing  the  values  of  /i  and  A  the  algorithm  was  tested  using  Narrow 
Band  Noise  signal  yh  for  its  performance.  This  signal  was  used  only  to  tune  the  value 
of  r  (i.e.  sharpness  of  the  notch  filter).  The  following  simulation  results  are  based  on 
parameter  adaptation  for  the  Dx(z)  given  by  equation  (4.17). 

5.  Response  for  Sinusoidal  signal 

A  sample  simulation  output  due  to  sinusoidal  signal  was  shown  in  Fig  [4.5 
-  4.7].  The  input  signal  is  composed  of  3  sinusoids  with  normalized  frequencies  ^/3, 
^/s,  and  H5/3  and  WGN  with  a  —  1.  The  spectrum  of  this  signal  is  shown  in  Fig 
4.5.  After  passing  this  signal  through  the  adaptive  filter  we  could  see  that  interferers 
were  removed  and  only  the  noise  was  left  behind.  This  is  clearly  shown  in  Fig  4.6. 
The  adaptation  process  is  shown  in  Fig  4.7. 

6.  Response  for  Narrow  Band  Noise 

In  this  testing  category  sinusoids  were  replaced  by  narrow-band  signals. 
These  signals  were  generated  via  equation  (4.4).  Superimposed  on  this  signal,  WGN 
with  cr=l  was  added.  A  typical  spectrum  of  this  signal  is  shown  in  Fig  4.8.  This 
signal  was  passed  through  the  adaptive  filter  and  these  narrow  band  interferers  are 
notched  out  and  only  WGN  is  left  out  as  shown  in  the  Fig  4.9.  Fig  4.10  shows  the 
corresponding  parameter  convergence. 
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7.  Tracking  FSK  Signal 

Transient  behavior  of  the  adaptive  filter  was  studied  by  applying  an  FSK 
signal  (generated  via  equation  4.11).  In  the  case  of  an  FSK  signal,  signal  spectrum 
constantly  changes.  This  property  is  useful  for  testing  the  tracking  ability  of  the 
adaptive  filter.  In  fact  a  WGN  with  a  =  1  was  also  added  to  the  signal.  An  adaptive 
filter  then  was  used  to  demodulate  the  signal,  by  tracking  the  spectrum.  This  is 
shown  in  Fig  4.13. 

8.  Suppressing  BPSK  Interference 

Having  seen  the  performance  of  the  adaptive  filter  under  various  conditions, 
it  is  only  needed  to  test  under  an  additional  constraint  of  broad  band  noise.  In  this 
case,  the  broad-band  signal  was  swamped  by  narrow-band  interferers  and  WGN.  A 
typical  broad-band  signal  is  generated  by  using  a  BPSK  signal  generator  via  equation 
(4.8).  In  this  we  have  used  carrier  at  0.25  cycle/sample  and  the  chip  frequency  at  0.125 
cycles/sample.  The  narrow-band  interference  signals  are  also  generated  by  using  the 
same  BPSK  signal  generator,  but  with  different  parameters.  Interferers  are  chosen 
such  that  the  chip  frequency  is  fixed  at  ^55/*  cycles/sample,  while  carriers  were  chosen 
a^  ifo-^51  36o^8'  an<^  iio^*  cycles/sample-  To  these  interferers  and  broad-band  signal, 
a  WGN  with  a  —  1  was  added.  The  composite  signal  spectrum  is  shown  in  Fig  4.14, 
indicating  the  three  interferers,  broad-band  signal,  and  WGN.  This  signal  was  passed 
through  the  adaptive  filter.  Output  of  the  filter  is  shown  in  Fig  4.15.  In  this  figure 
it  is  clear  that  the  interferers  were  removed  by  the  filter  and  only  the  desired  signal 
was  left  behind. 

9.  Model  order  mismatch 

The  model  order  of  the  algorithm  was  fixed  at  3  while  BPSK  signals  were 
generated  with  2  interferers  i.e.  m  =  2  and  Ar  =  3.  Fig  3.16  shows  the  parameter  con- 
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vergence  under  these  conditions  and  Fig  3.19  shows  for  m  =  1  and  N  =  3  conditions. 
We  could  notice  that  free  extra  parameter(s)  wandering  due  to  the  mismatch.  This 
has  the  effect  of  notching  out  the  desired  signal.  This  can  be  solved  by  deliberately 
introducing  a  known  BPSK  signal  at  fixed  place. 
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Summary  of  Tests  on  the  New  Algorithm 


TABLE  4.1 


Case 

sections 

Sines 

Vmin 

Umax 

Noise 

zero  order 

1st  order 

2nd  order 

V 

iter 

V 

iter 

V- 

iter 

la 

1 

1 

0.05 

20 

no 

0.05 

67 

0.05 

20 

0.05 

20 

Ian 

1 

1 

0.05 

20 

yes 

0.05 

jump 

0.05 

20 

0.05 

144 

lb 

1 

3 

0.05 

20 

no 

0.05 

jump 

0.05 

30 

0.05 

25 

lbn 

1 

3 

0.05 

20 

yes 

0.05 

jump 

0.05 

50 

0.05 

100 

2a 

5 

3 

0.10 

10 

no 

0.017 

340 

0.063 

200 

0.07 

400 

2an 

5 

3 

0.10 

10 

yes 

0.017 

1200 

0.063 

400 

0.07 

450 

2b 

5 

5 

0.10 

10 

no 

0.017 

1410 

0.063 

150 

0.07 

500 

2bn 

5 

5 

0.10 

10 

yes 

0.017 

3500 

0.063 

175 

0.07 

250 

2c 

5 

7 

0.10 

10 

no 

0.017 

1667 

0.063 

450 

0.07 

jump 

2cn 

5 

7 

0.10 

10 

yes 

0.017 

700 

0.063 

167 

0.07 

jump 

3a 

10 

7 

0.02 

50 

no 

0.02 

1667 

0.075 

500 

0.07 

1500 

3  an 

10 

7 

0.02 

50 

yes 

0.02 

5000 

0.075 

750 

0.07 

1000 

3b 

10 

10 

0.02 

50 

no 

0.02 

7500 

0.075 

2000 

0.07 

5250 

3bn 

10 

10 

0.02 

50 

yes 

0.02 

1000 

0.075 

2000 

0.07 

2500 

•  In  all  cases  with  noise,  parameters  converge  to  exact  value  when  number  of 
notches  is  greater  or  equal  to  number  of  sinusoids.  Iterations  listed  in  table 
represent  convergence  to  three  decimal  places. 

•  In  all  cases  with  noise,  parameter  oscillate  around  the  exact  value.  The  number 
of  iterations  listed  in  the  table  is  the  number  of  iterations  required  to  establish 
this  oscillatory  pattern. 

•  In  some  cases  with  more  interfering  sine  waves  than  notches  will  jump  at  random 
between  sine  waves.  This  pattern  is  indicated  in  the  table  by  the  word  jump. 

•  Values  of  vmin  and  vmax  in  the  table  refer  only  to  zero-order  forgetting  which 
generally  will  not  converge  without  limits  on  v.  All  other  types  of  forgetting 
were  run  without  limits  on  v 
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V.  HARDWARE  IMPLEMENTATION 

Although  hardware  design  is  beyond  the  scope  of  this  thesis,  in  this  chapter  we 
shall  briefly  look  at  possible  hardware  configurations  with  a  view  toward  feasibility 
of  the  new  algorithm. 

A.      HARDWARE  FEASIBILITY 

In  recent  years  many  dedicated  digital  signal  processing  chips  have  become 
commercially  available.  Table  5.1  gives  characteristics  of  some  typical  chips  that  are 
currently  available. 

TABLE  5.1 


Feature 

Commercial  Make 

MIPS  R3010 

Weitek  3364 

TI  8847 

AT&T  DSP32c' 

Cycle  Time(ns) 

40 

50 

30 

20 

Cycles/add 

2 

2 

2 

2 

Cycles/mult 

5 

2 

3 

2 

Cycles/divide 

19 

17 

11 

3 

Cycles/sq  root 

- 

30 

14 

? 

Data  on  MIPS  R3010,  Weitek  3364,  and  TI  8847  was  obtained  from  the  computer 
architecture  book  by  Hennesey  and  Patetrson  [Ref.  30]  while  the  data  for  AT&T 
DSP32c  is  obtained  from  DSP32c  data  manual  [Ref.  31].  Fig  5.1.  gives  a  schematic 
of  the  2nd  order  IIR  filter  hardware  scheme.  It  has  a  coefficient  memory  which  can 
be  set  by  an  external  device.  Close  examination  reveals  that  Multiplications  A  can 
be  performed  simultaneously  while  Multiplication  B  and  other  additions  are  to  be  in 
sequence.  Fig  5.2.  is  the  basic  building  block  for  the  adaptive  filter  and  is  identical  for 
all  the  sections.  This  is  offers  an  advantage  in  hardware  complexity  over  the  Kwan  k. 

Martin  approach.  Filters  1  and  2  have  transfer  functions  of  Hbp(z)  while  filter  3  has  a 

C<  WfeU^'d  do  Ja  • 
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transfer  function  Ha3(z).  Throughput  rate  is  primarily  determined  by  the  computing 
time  of  the  2nd  order  IIR  filter.  The  overall  block  diagram  of  the  adaptive  filter  for  a 
3-Notch  cancellation  was  given  in  Fig  5.3.  This  architecture  remains  same  either  for 
a  Floating  point  or  Fixed  Point  arithmetic. 
1.      Time  Budget 

In  the  architecture  of  Fig  5.1,  5.2,  and  5.3  the  most  vital  elements  are  the 
IIR  filter,  Controller  and  the  Forgetting  filter.  The  IIR  filter  corresponds  to  Hbp(z). 
The  controller  has  to  update  the  parameter  via  equation  (3.17)  which  calls  for  2 
multiplications,  1  addition,  and  1  division.  Similarly  the  Forgetting  Filter  has  to 
implement  equation  (3.19)  calling  for  2  multiplications  and  1  addition.  Results  are 
tabulated  as  given  below 

TABLE  5.2 


Element 

#  of  Effective 
Multiplications 

#  of  Effective 
Additions 

#of 

Divisions 

IIR  filter 

2 

3 

nil 

Controller 

2 

1 

1 

Fogetting 
Filter 

2 

1 

nil 

The  maximum  computing  time  is  at  the  IIR  filter.  Using  the  AT&T  DSP32c  proces- 
sor, it  takes  5x2  =  10  cylces  (ie  200ns)  for  implementation  the  desired  IIR  filter.  This 
corresponds  to  a  Throughput  rate  of  bMhz.  This  is  the  best  that  could  be  achieved 
with  the  existing  commercial  floating  point  processors.  However  processors  like  the 
HDSP66110  of  Honey  well  make  claim  an  ALU  speed  of  10ns  with  block  floating  point 
operations  can  conveniently  implement  at  \0Mhz  throughput  rate,  without  pipelin- 
ing the  IIR  filter.  With  pipelining  the  throughput  rate  would  increase  dramatically, 
but  will  never  exceed  the  maximum  throughput  rate  of  lOOMhz. 
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2.      VLSI  Approach 

Taking  advantage  of  the  fact  that  the  entire  hardware  can  be  generated  by 
a  simple  repetitive  use  of  the  building  block,  strong  consideration  should  be  given 
to  designing  a  VLSI  chip  for  Fig  5.2  rather  than  implementing  the  hardware  via  the 
existing  commercial  chips.  The  main  rationale  behind  this  idea  is  that  reliability 
of  the  system  is  very  important  in  EW  equipment.  Also  VLSI  has  the  potential  of 
ultimate-low  cost. 
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VI.  CONCLUSIONS 

In  this  thesis  the  problem  of  suppressing  narrow-band  interference  was  ad- 
dressed. The  problem  specific  to  the  Electronic  warfare  scenario  was  kept  in  mind 
while  solving  the  problem.  The  new  algorithm  [Ref.  26]  was  derived  and  simulated. 
This  new  algorithm  was  tested  against  various  signals  and  signal  conditions.  The 
results  are  highly  encouraging. 

Some  aspects  of  pipelining  were  discussed.  Reduced  hardware  complexity  of  the 
New  Algorithm  is  an  important  advantage  over  the  earlier  algorithm  by  Kwan  and 
Martin  [Ref.  8].  A  possible  hardware  scheme  for  implementing  the  new  algorithm 
was  discussed.  It  was  observed  that  with  the  existing  commercially  available  floating 
point  processors  a  throughput  rate  of  bMhz  is  achievable  while  using  processors  like 
HDSP66110  with  an  ALU  speed  of  10ns  it  is  possible  to  achieve  a  throughput  rate 
of  lOMhz. 

A.     FUTURE  WORK 

Future  directions  of  work  are 

i)  VLSI  design  of  the  system  with  an  architecture 

as  indicated  earlier  or  a  similar  architecture, 
ii)  Pipelining  of  2nd  Order  IIR  filter  suited  for 

this  application 
iii)  Design  using  recursive  DFT 

The  above  areas  of  work  are  the  logical  extensions  of  this  work.  However  a  radically 
different  approach  for  the  same  problem  using  adaptive  filtering  in  frequency  domain 
[Ref.  10]  should  also  be  considered. 
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APPENDIX    A 

dimension   faray(5000) 

dimension  uaray(5000) ,yaray(5000) 

dimension    fhat (5000) ,ph(5000) 

dimension   a(90),b(90) 

open  (unit=9 , f ile= ' spk. dat • , status= ' new • ) 
c 

c    This  programme  generates  input  output 
c    sequence  by  exciting  a  linear  system 
c    defined  by  the  numerator  polynomial 
c    B(z)  and  denominator  polynomial  A(z) 
c    this  data  is  stored  in  uarray  and  yarray 
c    Subsequently  an  ARMA  (pole, zero)  is  used 
c    to  fit  this  data, 
c 

n=1024 

ix=l 

yk=rand(ix) 

ix=0 

kk=8 
c 

c    generate  sequence  uk  and  yk 
c 

ip=3 

iz=2 

a(l)=1.7 

a(2)=-1.53 

a(3)=0.648 

b(l)=1.0 

b(2)=0.6 

do  10  k=l,n 

call  rbn(uk,ix,k) 

call  system(a,b,uk,yk, ip, iz ,k) 

uaray (k)=uk 

yaray (k)=yk 
10   continue 
c 

c    save  the  sequence  in  the  arrays  uaray  and  yaray 
c 

call   spktrm (yaray , faray) 
c 

c    spectrum  of  the  output  sequence  yaray  is  computed 
c    using  FFT  programme  and  stored  in  faray  for 
c    comparision 
c 

ip=4 

iz=l 
c 

c    Fit  an  ARMA(ip,iz)  model  for  the  given  data 
c    this  is  an  ordinary  least  squares  algorithm 
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print  *, * input-ip-iz — >' 

read  *,ip,iz 

do  2  0  k=l,n 

uk=uaray (k) 

yk=yaray (k) 

call   anna (a,b,uk,yk, ip, iz/kk/k) 
c  print   *, •num-> ■ , (b(i) , i=l, iz) 

c  print   *, 'den->' , (a(i) , i=l,ip) 

20        continue 

print   *, ' nui-> ' , (b(i) , 1=1, iz) 

print  *, 'den->' , (a(i) ,i=l, ip) 
c 
c    Sampling  of  the  Z-transform.  given 

c    coefficients  a  ....  and  b  

c    corresponding  H(z)  =  B(z)/A(z) 

c    is  evaluated  on  the  unit  circle 

c    for  obtaing  the  parametric  spectrum 


c 

c 
c 


call  zsmple (a,b, fhat ,ph, ip, iz) 


do  30  i=l,n/2 
c    print  *,i,' — > ' , fhat (i) , faray (i) 

write  (9,*)  i, fhat(i) ,faray(i) 
30   continue 

stop 

end 
c 
c 

subroutine  zsmple (a,b, r,ph, ip, iz) 

dimension  a(90),b(90) 

dimension  r (2048) ,ph (2048) 

complex  z ,ui, fn, fd, omega, delw, spec 

pi=atan(1.0) 

pi=4 .0*pi 

ui=(0. 0,1.0) 

delw=(pi/102  4.0)*ui*2.0 

omega=(0. 0,0.0) 

do  100  k=l,1024 

z=exp (-omega) 

fd=(0.0,0.0) 

do  10  i=l,ip 
10    fd=fd+a(i)*(z**(-i) ) 

fd=1.0-fd 

fn=(0. 0,0.0) 

do  20  i=l,iz 
20    fn=fn+b(i)*(z**(-i) ) 

spec=fn/fd 

r (k) =abs (spec) 

specr=real (spec) 
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speci=aimag (spec) 

ratio=speci/specr 

ph(k) =atand( ratio) 

omega=omega+delw 
c    print  *, ' — z-om-> ', z , omega 
100   continue 

return 

end 


c 
c 


c 
c 


function  atand(x) 

pi=atan(1.0) 

pi=4 . 0*pi 

xrad=atan(x) 

atand=xrad* (180. 0/pi) 

return 

end 


subroutine  rbn(b,ix,k) 
z=rand(ix) -0. 5 
if(z.ge.O.O)  b=1.0 
if(z.lt.O.O)  b=-1.0 
return 
end 

c 
c 

subroutine  system(a,b,uk,yk, ip, iz,k) 
dimension 

&    zkar (90) , zkma(90) , 

&     a(90),b(90) 

if  (k.ne.l)  go  to  250 

do  230  i=l,ip 
230   zkar(i)=0.0 

do  240  i=l,iz 
240   zkma(i)=0.0 
250   continue 

do  30  i=ip,2,-l 

30  zkar(i)=zkar(i-l) 
zkar(l)=yk 

do  31  i=iz,2,-l 

31  zkma(i)=zkma(i-l) 
zkma(l)=uk 

c 
c 

yk=0 . 0 

do  200  i=l,ip 
200  yk=yk+a(i) *zkar(i) 

do  210  i=l,iz 
210  yk=yk+zkma(i) *b(i) 
return 
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end 


c 
c 


subroutine  delay  (uk,yk,  idelay,k) 

dimension  d(500) 

if  (k.ne.l)  goto  10 

do  20  i=l,500 
20   d(i)=0.0 
10   continue 

do  30  i=500,2,-l 
30   d(i)=d(i-l) 

d(l)=uk 

yk=d(idelay+l) 

return 

end 
c 
c 

subroutine  spktrm(taray, faray) 

complex  u(5000) 

dimension  taray(5000) , faray (5000) 

m=10 

n=1024 

rn=n 

do  20  i=l,n 

u(i)=taray (i) 
20    continue 

call  pfft(u,m,n) 

do  30  i=l,n/2 

faray (i)=real (u(i) ) 
30   continue 

return 

end 
c 
c 
c    ********************************** 

c  * 

c  Subroutine  for  computing  DFT  of     * 

c  an  array  'a'  is  complex  and  a  pair 

c  of  numbers  are  to  be  specified      * 

c  for  each  point  * 

c  m  is  the  2  power  index         * 

c  say  m=10  then  number  of  points   * 

c  are  1024  * 

c  after  computation  fft  is  kept    * 

c  in  the  same  complex  array  'a'     * 

c  * 

c  ********************************** 

c 

subroutine  pfft(a,m,n) 
complex  a (5000) ,u,w, t 
complex  ui,ur 
pi=3. 14 1592 653 589793 
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n=2**m 

nv2=n/2 

nml=n-l 

j=l 

do  10  i=l,nml 

if(i.ge.j)  goto  20 

t=a(j) 

a(j)=a(i) 

a(i)=t 
20    k=nv2 
30    if(k.ge.j)  goto  10 

j=j-k 

k=k/2 

goto  3  0 
10   j=j+k 

do  40  1=1, m 

le=2**l 

lel=le/2 

u=(l. 0,0.0) 

w=cmplx(cos (pi/lel) , -sin(pi/lel) ) 

do  40  j=l,lel 

do  50  i=j,n,le 

ip=i+lel 

t=a (ip) *u 

a(ip)=a(i)-t 
50   a(i)=a(i)+t 
4  0   u=u*w 

rn=n/2 

ui=(0.0, 1.0) 

ur=(1.0, 0.0) 

do  60  i=l,n/2 

a(i)=a(i)/rn 

a  ( i+512 )=a( i+512 )/rn 

areal=real (a(i) ) 

aimgn=aimag(a (i) ) 

amagn=abs ( a ( i ) ) 

aphase=atand(aimgn/areal) 

a (i)=ur*amagn+ui*aphase 

areal=real (a (i+512) ) 

aimgn=aimag(a (i+512) ) 

amagn=abs(a(i+512) ) 

aphase=atand(aimgn/areal) 

a ( i+512 ) =ur*amagn+ui*aphase 
60   continue 

return 

end 
c 
c 

c    $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 
c 
c 

subroutine  arma (a,b,uk,yk, ip, iz,kk,k) 
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dimension 

&    zkar (90) , zkma(90)  , 

&     zk(90) ,a(90) ,b(90) ,c(90) 

if  (k.ne.l)  go  to  250 

n=ip+iz 

ipp=ip+l 

do  230  i=l,ip 
230   a(i)=0.1 

do  240  i=l,iz 
240  b(i)=0.1 
250   continue 

do  30  i=ip,2,-l 

30  zkar (i)=zkar (i-1) 
zkar (l)=ykold 
ykold=yk 

do  31  i=iz,2,-l 

31  zkma (i)=zkma (i-1) 
zkma (l)=uk 

do  32  i=l,ip 
c(i)=a(i) 

32  zk(i)=zkar (i) 
do  33  i=l,iz 
c(i+ip)=b(i) 

33  zk(i+ip)=zkma(i) 

call  syseq(c,yk, zk, n, kk,k) 
do  34  1=1,1? 

34  a(i)=c(i) 

do  35  i=l,iz 

35  b(i)=c(i+ip) 

return 
end 
c 
c 
c 

subroutine  syseq(a,yk, zk,n,kk,k) 
c 
c    *************************************** 

c  solves  system  of  equations  with 

c  #  of  equations  more  than  unknowns 

c  using  linear  reggression. 

c  t 

c  equation  is  yk=zk  *  a 

c  where  'a'  is  n  vector  to  be  estimated. 

c 

c  **************************************** 

c 

dimension 
&     del(90) ,phat(90,90) ,zk(90) ,y(90) ,a(90) 
if  (k.ne.l)  go  to  250 
alpna=l. 0 
250   continue 

call  inv(phat, zk, alpha, n,kk,k) 
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ykhat=0. 0 
do  40  i=l,n 
40   ykhat=ykhat+zk(i)*a(i) 
er=yk-ykhat 

do  60  i=l,n 
60   y(i)=er*zk(i) 

do  80  i=l,n 

del(i)=0.0 

do  80  j-l,n 
80   del(i)=phat(i,j)*y(j)+del(i) 

if (mod(k,kk) .ne.O)  goto  70 

c    print  *, • >' ,  k 

c    print  100,  ( (phat (i, j ) , i=l,n) , j=l,n) 

c    print  130,  (del (i) , i=l,n) 

c    print  160, er 

c    print  170,  (zk( j) , j=l,n) 

70   continue 

do  50  i=l,n 
50   a(i)=a(i)+del(i) 
100   format (2x, 'phat' ,2x,/,4el6.  6) 
120   format(2x, 'zk' , 2x,/,4el6. 6) 
130   format(2x, 'del' ,2x,/,4el6.6) 
160   format(2x, 'ekl' ,2x,el6.6,/) 
170   format(2x, 'zk' ,2x,/,4el6.6) 


c 
c 


return 
end 


subroutine  inv(phat, zk, alpha, n,kk,k) 

dimension  phat(90,90) ,zk(90) ,q(90) ,r(90) ,delp(90,90) 

if(k.ne.l)  goto  5 

do  6  i=l,n 

do  6  j=l,n 

if(i.eq.j)  phat (i, j ) =1. 0 

if(i.ne.j)  phat (i, j )=0. 0 
6    continue 
5    continue 

do  90  i=l,n 

do  90  j=l,n 

phat ( i , j ) =phat ( i , j ) /alpha 
90   continue 

do  10  i=l,n 

q(i)=0.0 

do  10  j=l,n 
10   q(i)=phat(i,j)*zk(j)+q(i) 
sum=1.0 

do  20  i=l,n 
2  0   sum=sum+zk(i) *q(i) 

do  30  j=l,n 

r(j)=0.0 

do  30  i=l,n 
30   r(j)=phat(i,j)*zk(i)+r(j) 
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do  40  j=l,n 

do  40  i=l,n 
40   delp(i,j)=q(i)*r(j) 

do  50  j=l,n 

do  50  i=l,n 
50   phat ( i ,  j  ) =phat ( i , j ) - (delp ( i , j ) /sum) 

return 

end 
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APPENDIX  B 
c 

c    simulation  of  michal  paper  IEEE 
c 

dimension  akl (10) ,ak2 (10) ,s(10) ,xout(10) , theta (10) 

dimension  aks (10) ,wp(10) ,ak3 (10) 

dimension  array (4 , 4000) 

open (  unit=9 , f ile= ' martin . dat ' , status= • new ' ) 

open (  unit=8 , f ile= • mart . dat ' , status= ' new ' ) 

print  *, • -input-SNR-in-dB >' 

read  *,snr 
c    if  (snr.ne.0)  goto  200 
c    var=0.5 
c    goto  210 

var=1.0/(2.0*(10.0**(snr/10.0) )) 
210   continue 

print  *,' — variance — >',var 

n=l 

print  *,'-input-nh — >• 

read  *,n 

print  *, * -input-#-of-waves — >' 

read  *,nd 

print  *, '-input-angles — >' 

read  *, (theta (i) , i=l,n) 

r=0.9 

pi=4 .0*atan(1.0) 

rad=pi/180.0 

amu=0. 0 
c    do  50  kl=l,50 

amu=0.015 

print  *,' -input-step-size —  0.015 — >',amu 

read  * , amu 

print  *, • — angles — > ' , (theta (i) , i=l,n) 

print  *, • — input — 1-for-step-change — else-0 — >' 

read  *,ic 

if  (ic.eq.l)  print  *,' — input-step-in-degrees — >' 

if  (ic.eq.l)  read  *,step 

do  40  i=l,n 

ak2(i)=(l-r*r) 

ux=theta ( i ) *rad 

value=cos (ux) 
c    akl(i)=dl*sqrt(1.0+r*r-2.0*r*value) 

akl ( i ) =sqrt ( 1 . 0+r*r ) 

aks ( i ) =sqrt ( 1 . 0+r*r-2 . 0*r*value) 
40   continue 
c    amu=amu+ 0.0001 

do  10  k=l,6500 

call  system ( yk ,wk, theta, step ,var, ic,ndf k) 

call  f ilters(akl,ak2 , aks, s,en,xout,ykf n, k) 

cal 1  increment ( akl , ak2 , aks , s , en , amu , avg , n , k) 

do  30  i=l,n 
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xx=2.0*sqrt(l-0.5*ak2(i) ) 

wp(i)=2.0*asin(akl(i)/xx) 
c    wp(i)=wp(i)/rad 

wp(i)=cos(wp(i) ) 
3  0   continue 

do  60  i=l,n 

ak3(i)=abs(akl(i) ) 
60   continue 

kk=mod(k,100) 

ll=mod(k,10) 

if(kk.eq.O)  print  *, • — akl — >' , (ak3 ( j ) , j=l,n) 

if(ll.eq.O)  write(9,*)  k, (ak3 ( j ) , j=l,n) 

if(kk.eq.O)  print  *, ' — aks — >' , (aks ( j ) , j=l,n) 
c    if(kk.eq.O)  print  *,k, '-avg->' ,  avg 
c    print  *, (wp(i) ,i=l,n) 
c    write^,*)  k,  (wp(i)  ,  i=l,n) 

if (k.gt. 1500. and. k. It. 2500)  writefS,*)  k,yk,en 

if (k.gt. 5000. and. k. It. 6000)  write(8,*)  k,yk,en 
10   continue 

c  print  *,  kl,amu,avg 
c  write  (9,*)  kl^amu^vg 
50   continue 

close (9) 

stop 

end 


subroutine  increment ( akl , ak2 , aks ,s,en, amu , avg , n , k) 

dimension  akl (10) ,ak2 (10) , aks (10) ,s(10) ,ss(10) 

dimension  dec(10) ,flag(10) ,fading(10) 

if  (k.ne.l)  goto  30 
c    amu=0.001 

avg=0. 0 

fading(l)=0.9 

fading(2)=0.9 

fading(3)=0.9 

fading(4)=0.9 
30   continue 

do  10  i=l,n 

ss(i)=fading(i)*ss(i)+s(i)*s(i)*(1.0-fading(i) ) 
10   continue 
c    do  31  i=l,n 

c31   print  *, '-en-s-ss->' ,en,s(i) ,55(1) 
c    sum=0.0 

do  20  i=l,n 

dec (i)=amu*en*s(i)/(ss(i) +0.0001) 

akl(i)=akl(i)-dec(i) 
c    sum=sum+ ( abs ( aks ( i ) -akl ( i ) ) ) / ( aks ( i ) ) 

call  stability (akl, ak2, flag,n) 

akl(i)=akl(i)+flag(i)*dec(i) 
2  0   continue 

realk=k 
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c 
c 
c 


c 
c 


avg=avg* ( (realk-1. 0)/realk) +sum/realk 

return 

end 


subroutine  filters (akl, ak2 , aks, s, en,xout,xin,n,ki) 

dimension  w(10,3) ,a(10,2) ,b(10,2) ,6(11) ,xbp(10) 

dimension  g(10, 3) ,gbp(10) ,u(10) 

dimension  q(10,3) ,s(10) ,xout(10) ,c(10) 

dimension  akl (10) ,  ak2 (10) , aks (10) ,35(10) 
c 

if  (ki.ne.l)  goto  200 
c    open(  unit=9 , file=' martin. dat1 , status=' new' ) 
200   continue 

i=n 

do  40  i=l,n 

a(i,l)=-(akl(i)*akl(i)+ak2(i)-2.0) 

as(i)=-(aks(i)*aks(i)+ak2(i)-2.0) 

a(i,2)=-(1.0-ak2(i) ) 

b(i,l)=-0.5*ak2(i) 

b(i,2)=0.5*ak2(i) 

c(i)=-akl(i)*ak2(i) 
4  0   continue 
c    do  41  i=l,n 

c    print  *, '-kl-k2->' ,akl(i) ,ak2(i) 

c4l  print  *, '-a-b-c-^ ,a(i, 1) ,a(i,2) ,b(i, 1) ,b(i,2) ,c(i) 
c 
c 

e(l)=xin 

i=n 

do  10  j=l,n 

w(i,3)=w(i,2) 

w  ( i ,  2 ) =w ( i  1 ) 

w(i,l)=a(i,l)*w(i,2)+a(i,2)*w(i,3)+e(j) 

xbp(j)=b(i,l)*w(i,l)+b(i,2)*w(i,3) 

e(j+l)=xbp(j)+e(j) 

i=i-l 
10   continue 

en=e (n+1) 


i=n 
do  20  k=l,n 
u(k)=en-xbp(k) 

g(i,3)=g(i,2) 

g(i,2)=g(i,l) 

g(i,l)=a(i,l)*g(i,2)+a(i,2)*g(i,3)+u(k) 

gbp(k)=b(i,l)*g(i,l)+b(i,2)*g(i,3) 
q(k,3)=q(k,2) 

q(k,2)=q(k,l) 
q(k,l)=a(i,l)*q(k,2)+a(i,2)*q(k,3)+gbp(k) 

79 


s(i)=c(i)*q(k,2) 
i=i-l 
20   continue 

kk=mod(ki,  10) 
c    if(kk.eq.O)  write(9,*)  ki,en,a(l,l) ,as(l) , a(2,l) , as(2) 
c    if(kk.eq.O)  print  *,    ki, (a (i, 1) , i=l,n) 
c    if(kk.eq.O)  print  *,    ki, (as (i) ,  i=l,n) 

c    if (ki.ge.1800)  write(9,*)  ki,e(l) ,en,xbp(l) ,3(1, 1) ,as (1) 
c    if (ki.ge. 1800)  print  *,    ki,e(l) ,en,xbp(l) ,a(l, 1) ,as (1) 

return 

end 
c 
c 
c 

subroutine  system  (yk,wk,  theta,  step,var,  ic,n,k) 

dimension  theta(lO) ,yout(10) ,yz(10) 

if  (k.ne.l)  goto  100 
100  continue 

c    call  bpf (theta, yz ,n,k) 
c    call  waves  (theta,yout,n,k) 

call  f ile(yk,theta,step, ic,n,k) 

sum=0. 0 

wk=0.0 

do  10  i=l,n 

sum=sum+yout ( i ) 

wk=wk+yz (i) 
10   continue 

call  gnoise (g, var , k) 
c    yk=sum+g 
c    yk=wk+g 

yk=yk+g 

return 

end 
c 
c 

subroutine  waves (theta, yout,n,k) 

dimension  zkar(10,3) ,theta(10) ,yout(10) ,al(10) 

if(k.ne.l)  goto  10 

pi=(atan(1.0) )*4.0 

rad=pi/180. 0 

do  20  i=l,n 

zkar(i, 1)=0.0 

zkar (i,2)=-sin(theta(i) *rad) 

al (i) =2 . 0*cos (theta ( i) *rad) 

2  0   continue 
10   continue 

do  30  i=l,n 

zkar ( i , 3 ) =zkar ( i ,  2 ) 

zkar (i, 2)=zkar (i,  1) 

zkar(i, l)=al(i) *zkar (i,2) -zkar (i, 3) 

yout (i)=zkar(i,l) 

3  0   continue 
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return 
end 


c 
c 


subroutine  gnoise (g, s,k) 

if(k.ne.l)  goto  10 

ix=l 

yf l=rand(ix) 

ix=0 
10   continue 

sum=0. 0 

do  20  i=l,12 

yf l=rand(ix) 
2  0   sum=sum+yfl 

g= (sum-6. 0) *s 

return 

end 
c 
c 

subroutine  stability (akl, ak2 , flag,n) 

dimension  akl (10) ,ak2 (10)  , flag (10) , a(10,2) 

do  30  1*1,  n 

flag(i)=0.0 

a(i,l)=-(akl(i)*akl(i)+ak2(i)-2.0) 

a(i,2)=-(1.0-ak2(i) ) 

des=a(i,l) *a  (i,  1) +4 . 0*a (i, 2) 

if  (des.lt. 0.0)  goto  10 

des=0. 5*sqrt (des) 

rootl=0.5*a(i, l)+des 

rl=abs (rootl) 

root2=0.5*a(i, 1) -des 

r2=abs(root2) 

if  (rl.ge.1.0)  goto  20 

if  (r2.ge.l.O)  goto  20 
c    print  *, ' -real-roots — >•, rootl, root2 

goto  30 
20   flag(i)=1.0 

goto  3  0 
10   des=-des 

y=0. 5*sqrt (des) 

x=0.5*a(i,l) 

xx=abs(x) 

if  (xx.le. 1.0e-6)  x=0. 0000001 

angle=(atan(y/x) ) *57 . 3 

radius=sqrt (x*x+y*y ) 
c    print  *, ' — complex — roots — > ' , radius, angle 

if (radius. ge. 1. 0)  goto  20 
30   continue 

return 

end 
c 
c 
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$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 


c 
c 


subroutine  bpf( theta, yout,n,k) 

dimension  zkar(10,3)  ,a(10,3)  ,zkma(10,4)  ,b(10,3) 

dimension  theta(lO) ,g(10) ,yin(10) ,yout(10) ,uk(10) 

if  (k.ne.l)  go  to  10 

pi=4.0*atan(1.0) 

rad=pi/180.0 

do  15  i=l,n 

zkar (i, 1)=0. 0 
15    zkar(i,2)=0.0 

r=0.985 
c    print  *,  ' — input — r — >',r 
c    read  *,r 

do  25  i=l,n 

a(i/l)=l 

a(i,3)=r*r 

g(i)=0.5*(l-r*r) 

b(i,l)=g(i) 

b(i,3)=-g(i) 

a ( i , 2 ) =2  *r*cos (theta ( i) *rad) 
25    b(i,2)=0.0 
10   continue 

save=var 

do  35  i=l,n 

var=2 . 0 

call  gnoise (gk, var,k) 

yin(i)=gk 

zkar ( i , 3 ) =zkar ( i ,  2 ) 

zkar (i, 2) =zkar (i,  1) 

zkma ( i ,  3 ) =zkma ( i , 2 ) 

zkma (i,2)=zkma(i,  1) 

zkma (i, l)=yin(i) 

uk ( i ) =zkma ( i , 1 ) *b ( i , 1 ) -zkma ( i , 2 ) *b ( i , 2 ) +zkma ( i , 3 ) *b ( i , 3 ) 

zkar(i, l)=zkar(i,2) *a(i, 2) -zkar (i, 3) *a(i, 3)+uk(i) 
c    print  *, • — fil — >• ,afil, theta ,r 
c    print  *,' — fil — >',zkar 
35   yout(i)=zkar (i,  1) 

var=save 
return 
end 
c    $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 
c 

subroutine  f ile (yk, theta, step , ic,n,k) 

dimension  output (8000) , theta (10) 

if  (k.ne.l)  goto  10 

call  gendata (output, theta, step, ic,n,k) 
10   continue 

yk=output (k) 

return 

end 
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c  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
c  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
c      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


subroutine  gendata (output, theta, step, ic,nn,k) 

DIMENSION  YSAMP(65536) , data (4 , 8000) , output (8000) 

dimension  theta (10) 
C     OPEN  (UNIT=9,FILE='bpsk.dat' ,STATUS=,NEW) 

n=7000 

amag=1.414 

tcarr=4 . 8 

tchip=7.0 

ichip=0 

tdata=n 

tdelay=0. 0 

1=9999 

tchip=200.0 

do  10  i=l,nn 

tcarr=3  6  0 . 0/theta ( i ) 

call  bpsk(n,amag, tcarr, tchip, ichip, tdata,tdelay, l,ysamp) 

do  20  j=l,n 

if  (ic.eq.0)  goto  80 

if  (i.eq.2)  goto  50 

data ( i , j ) =y samp ( j ) 

goto  60 
50    if  (j.gt.3000)  ysamp(j)=0.0 
80   data(i, j )=ysamp( j ) 
60   continue 
20   continue 

print  *,i,' — > ' ,theta (i) ,tcarr 
10   continue 

tcarr=360. 0/90.0 

tchip=8.0 

call  bpsk(n,amag, tcarr, tchip, ichip, tdata,tdelay, l,ysamp) 

do  30  i=l,n 

x=0.0 

do  40  j=l,nn 
4  0   x=x+data(j , i) 

output ( i ) =x+ysamp ( i ) 
c    write (9,*)  x 
30   continue 

if  (ic.eq.0)  return 

nl=n-3000 

tchip=200.0 

thetal=theta ( 2 ) +step 

tcarr=3  60 . 0/thetal 

cal 1  bpsk (nl , amag , tcarr , tchip , ichip , tdata , tdelay , 1 , ysamp ) 

j=l 

do  70  i=3000,n 
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output ( i ) =output ( i ) +ysamp ( j ) 

j=j  +  l 
7  0   continue 
c    close(9) 

return 

end 


this  programme  was  modified  to  suit  the  martin 

programme 

subroutine  bpsk 

&    (n , amag , tcarr , tchip ,  ichip , tdata , tdelay ,11, ysamp) 
c 
c 
c 

C     THIS  PROGRAM  GENERATES  SAMPLES  OF  A  DIRECT  SEQUENCE  BI-PHASE 
SHIFT 

C     KEYED  SPREAD  SPECTRUM  SIGNAL.   THE  "INFORMATION"  BITS  AND  THE 
C     SPREADING  SEQUENCE  USED  ARE  RANDOMLY  GENERATED.   PARAMETERS 
REQUIRED 


C 
C 
C 

C 
C 
C 
C 
C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

THE 
C 


FOR  OPERATION  ARE: 

N   =   NUMBER  OF  SAMPLES  GENERATED 

FOR  CONSISTENCY  WITH  USE  BY  AN  FFT 
ALGORITHM,  N  SHOULD  BE  AN  INTEGER 
POWER  OF  2  —  TYPICALLY  1024 

NOTE:  IF  N>1024  DIM  OF  YSAMP  MUST  CHANGE 
MAG      =   MAGNITUDE  OF  CARRIER  WAVEFORM 
TSDELAY  =   NUMBER  OF  SAMPLE  TIMES  DELAYED  FROM 
t  =  0  BEFORE  BEGINNING  SAMPLES 
AN  ARBITRARY  VALUE  THAT  ALLOWS  SOME 
FLEXIBILITY  IN  SAMPLING.   SHOULD  BE 
ZERO  FOR  SAMPLING  AT  t=0. 
TDATA    =   NUMBER  OF  SAMPLE  TIMES  IN  ONE  DATA  BIT 
TCHIP    =   NUMBER  OF  SAMPLE  TIMES  IN  ONE  BIT  OF 

SPREADING  CODE 
TCARR    =   NUMBER  OF  SAMPLE  TIMES  IN  ONE  CYCLE  OF 

CARRIER  FREQUENCY 
TSAMP    =   DURATION  OF  A  SAMPLE  TIME.   IN  GENERAL 
TSAMP  WILL  ALWAYS  BE  =  1.0,  SINCE 
TIME  SAMPLED  VALUES  BECOME  STATIC 
WHEN  STORED  AND  CAN  BE  SCALED  LATER. 

THE  ALGORITHM  USED  TO  GENERATE  THE  DS-BPSK  SIGNAL  RELIES  UPON 

BUILTIN  FORTRAN  RANDOM  NUMBER  GENERATOR  RAND(L)  TO  PRODUCE  THE 


C     INFORMATION  BIT  STREAM  AND  THE  SPREADING  SEQUENCE  CODE.   AT 

THE  TIME 

C     OF  EACH  SAMPLE,  THESE  TWO  BITS  ARE  MODULO  2  ADDED  TO  PRODUCE 

A  CHIP 

C     BIT.   THIS  BIT  IS  THEN  USED  TO  SHIFT  THE  PHASE  OF  THE  CARRIER 

BY 
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C     18  0  DEGREES  IF  1  OR  BY  0  DEGREES  IF  0.   THE  RESULTANT  VALUE  OF 

THE 

C     COSINE   FUNCTION   IS   MULTIPLIED   BY   THE   AMPLITUDE   OF   THE 

WAVEFORM,  THEN 

C     STORED  IN  FILE  BT.DAT  IN  THIS  FILE  DIRECTORY.   THE  PROCESS  IS 

CONTIN- 

C     UED  UNTIL  N  SAMPLES  ARE  GENERATED.   NOTE  THAT  FIRST  LINE  OF 

BT . DAT 

C     CONTAINS  N,  DSBPSK,  DATE  AND  TIME  AS  ID  FIELD  FOR  FILE. 

C 

C  W.R.TUCKER  9-MAY-83 

C  CONVERTED  PROGRAM  TO  ALLOW  FOR  M-SEQUENCE  CHIP 

C  M-SEQUENCE  GENERATOR  REQUIRES  DATA  FILE  FSR.DAT  WITH 

C  PARAMETERS  FOR  FEEDBACK  SHIFT  REGISTER. 

C  W.R.TUCKER  4-AUG-83 

C 

C     INITIALIZE  —  I  =  INFO  BIT  NUMBER,  J  =  SP  SEQ  CODE  BIT  NUMBER 

C     L  IS  ARGUMENT  FOR  RAND(L) — A  DIFFERENT  SEQUENCE  MAY  BE 

C     GENERATED  BY  INITIALIZING  L  WITH  DIFFERENT  VALUES. 

C 

C    Modified  by  HHL  for  installation  on  ULTRIX  ECE  VAX,  12/90 

c    Writes  output  in  file  bpsk.dat 

c 

DIMENSION  YSAMP(65536) 

INTEGER  I,J,L, INFO, CHIP, IFSR, IFBD (100) ,IVAL(100) 

REAL  TDATA,TCHIP,TCARR,TDELAY,TSAMP,MAG,YARG,YSAMP,PI 

PI  =  3.14159265 

TSAMP  =1.0 

1  =  0 

J  =  0 

L  =  0 

WRITE(6,900) 
900   FORMATC   GENERATION  OF  DIRECT  SEQUENCE  BPSK  SIGNAL', 

A  •  IN  FILE  bpsk.dat') 

WRITE(6,1000) 

1000  FORMATC   #  SAMPLES  TO  GENERATE        =  ',$) 
c     READ ( 5 ,  * )  N 

print  *,n 
WRITE(6,1001) 

1001  FORMATC   MAX  AMPLITUDE  OF  SAMPLE  VALUES  (R)=  ',$) 
C     READ(5,  *)  MAG 

mag=amag 
print  *,mag 
WRITE(6,1002) 

1002  FORMATC   #  SAMPLES  PER  CARRIER  CYCLE  (R)=  ',$) 
C     READ (5,  *)  TCARR 

print  *,tcarr 
WRITE (6, 1003) 

1003  FORMATC   #  SAMPLES  PER  CHIP  BIT    (R)  =  ',$) 
C     READ (5,  *)  TCHIP 

print  *,tchip 
WRITE (6, 102  0) 
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1020  FORMAT (•    ENTER  (0):RANDOM  CHIP,  OR  (1):REPEAT  M-SEQ  CHIP 

=  l  t  $) 

C     READ(5,  *)  ICHIP 

print  *,ichip 

WRITE (6, 1004) 

1004  FORMAT ('   #  SAMPLES  PER  INFO  BIT    (R)  =  ',$) 
C     READ (5,  *)  TDATA 

print  *,tdata 
WRITE(6,1005) 

1005  FORMAT ('   #  DELAYS  BEFORE  SAMPLING    (R)=  ',$) 
c     READ (5,  *)TDELAY 

print  *,tdelay 
WRITE(6,1006) 

1006  FORMATC   RANDOM  NUMBER  SEED  (14) >  ',$) 

1=11 

C     READ (5,  *)L 

print  *,1 
c    Initialize  RAND  by  calling  with  input  non-zero  L. 
c    Subsequent  calls  will  be  with  L  =  0. 

RANDOM=RAND(L) 
c    debug 
C     WRITE (6,*)  '  SEED  AND  RANDOM  NUMBER  RETURNED: '  , L, RANDOM 

L=0 

WRITE(6,1010) 
1010  F0RMAT(/1X,  ' WORKING ') 

IF  (ICHIP  .EQ.  0)  GO  TO  5 

OPEN  (UNIT  =  1,  FILE  =  ' FSR.DAT',  STATUS  =  'OLD') 

READ  (1,500)  IFSR 

READ  (1,600)  (IFBD(K)  ,  K  =  1,  IFSR) 
500   FORMAT (13) 
600   FORMAT(I3) 

CLOSE  (UNIT  =  1) 

DO   5  K  =  1  ,  IFSR 

IVAL(K)  =  1 
5     CONTINUE 

DO  100  K  =  1,N 
C   CHECK  TO  SEE  IF  WE  NEED  TO  GENERATE  A  NEW  DATA  BIT 

IF  ( (K+TDELAY) *TSAMP  .LT.  I*TDATA*TSAMP)  GO  TO  10 
C   IF  SO  DO  IT  HERE 

RANDOM  =  RAND(L) 
c    debug 
C     WRITE (6,*)  '  IFLAG  AND  RANDOM  NUMBER  RETURNED: •  , L, RANDOM 

INFO  =  0 

IF  (RANDOM  .GT.  0.5)  INFO  =  1 
C   KEEP  TRACK  OF  WHICH  DATA  BIT  WE  ARE  ON 

I  =  1+1 
10    CONTINUE 
C   NOW  CHECK  TO  SEE  IF  WE  NEED  TO  GENERATE  A  NEW  CHIP  BIT 

IF  ( (K+TDELAY) *TSAMP  .LT.  J*TCHIP*TSAMP)  GO  TO  20 

IF  (ICHIP  .NE.  1)  GO  TO  15 

CHIP  =  IVAL(IFSR) 

CALL  MSEQ(IFSR,IFBD,IVAL) 
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GO  TO  19 
15  CONTINUE 
C   IF  SO  DO  IT  HERE 

RANDOM  =  RAND(L) 

CHIP  =  0 

IF  (RANDOM  .GT.  0.5)  CHIP  =  1 
C   KEEP  TRACK  OF  THE  CHIP  BIT  NUMBER 
19    CONTINUE 

J  =  J+l 
2  0    CONTINUE 
C   NOW  WE  DETERMINE  THE  PHASE  SHIFT 

PHASE  =  FLOAT (MOD (INFO  +  CHIP, 2)) 
C  COMPUTE  THE  ARGUMENT  FOR  THE  COSINE  FUNCTION 

YARG  =  2.0*PI*(1.0/(TCARR*TSAMP) )*(K+TDELAY)   +  PI*PHASE 

YSAMP(K)=  MAG  *  COS (YARG) 
100   CONTINUE 
C  NOW  SAVE  THE  VALUE 

c   ranga   OPEN  (UNIT=1, FILE='bpsk.dat ' ,STATUS= 'NEW' ) 
C     CALL  DATE(BDATE) 
C     CALL  TIME(CTIME) 

c-ranga    WRITE ( 1 , 12  5 ) N , MAG , TCARR , TCHIP , TDATA 
125      FORMAT ( I 5 , 5X , F4 . 1 , 10H*DSBPSK+M2 , F7 . 2 , lHf , F7 . 2 , lHc , F7 . 1, 

*    lHd) 
c-ranga    WRITE (1,150)  (YSAMP(I) , 1=1, N) 

do  121  i— 1,11 

c    write (1,*)  ysamp(i) 
121   continue 

150  FORMAT (8F16. 8) 
CLOSE  (UNIT=1) 
return 

END 

SUBROUTINE  MSEQ (IFSR, IFDB, IVAL) 
C     THIS  SUBROUTINE  PERFORMS  THE  SHIFTING  OPERATION  OF  AN  IFSR 
STAGE 

C     FEEDBACK   SHIFT   REGISTER,   WITH   FEED   BACK   CONNECTIONS   AS 
INDICATED 

C     BY  IFDB.    IVAL  IS  THE  INITIAL  CONTENTS  OF  THE  FSR  AND  WILL 
CONTAIN 

C     THE  FINAL  CONTENTS  AFTER  SHIFTING. 
C  W.R.  TUCKER  4  AUG  83 

INTEGER  IFDB(IFSR) ,  IVAL(IFSR) 

ISUM  =  0 

DO  10  I  =  1,  IFSR 

ISUM  =  IFDB(I)  *  IVAL(I)  +  ISUM 
10    CONTINUE 

IBIT  =  MOD(ISUM,2) 

DO  20  I  =  1,  IFSR  -  1 

IVAL (IFSR  +1-1)  =  IVAL (IFSR   -  I  ) 
2  0    CONTINUE 
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IVAL(l)  =  IBIT 

RETURN 

END 

C  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

c  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

C  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

C  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
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